multiple stressor
Marcel, Owen
12/2/2021
Introduction
This markdown contains the full-length 16S rRNA amplicon sequencing analysis of 240 micro-ecosystems, which were treated with 4 stressors in all possible combinations at three different temperatures. The analysis is related to the manuscript “Stratified microbial communities are highly sensitive towards multiple combined global change factors, revealing antagonistic and synergistic effects” by Suleiman et al.
library(Biostrings)## Loading required package: BiocGenerics## Loading required package: parallel##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which, which.max, which.min## Loading required package: S4Vectors## Loading required package: stats4##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:base':
##
## expand.grid## Loading required package: IRanges## Loading required package: XVector##
## Attaching package: 'Biostrings'## The following object is masked from 'package:base':
##
## strsplitlibrary(ShortRead)## Loading required package: BiocParallel## Loading required package: Rsamtools## Loading required package: GenomeInfoDb## Loading required package: GenomicRanges## Loading required package: GenomicAlignments## Loading required package: SummarizedExperiment## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.## Loading required package: DelayedArray## Loading required package: matrixStats##
## Attaching package: 'matrixStats'## The following objects are masked from 'package:Biobase':
##
## anyMissing, rowMedians##
## Attaching package: 'DelayedArray'## The following objects are masked from 'package:matrixStats':
##
## colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges## The following objects are masked from 'package:base':
##
## aperm, apply, rowsumlibrary(reshape2)
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:Biobase':
##
## combine## The following object is masked from 'package:BiocGenerics':
##
## combinelibrary(phyloseq)##
## Attaching package: 'phyloseq'## The following object is masked from 'package:SummarizedExperiment':
##
## distance## The following object is masked from 'package:Biobase':
##
## sampleNames## The following object is masked from 'package:GenomicRanges':
##
## distance## The following object is masked from 'package:IRanges':
##
## distancelibrary(boot)
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.3 ✓ purrr 0.3.4
## ✓ tibble 3.1.0 ✓ dplyr 1.0.4
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::collapse() masks Biostrings::collapse(), IRanges::collapse()
## x dplyr::combine() masks gridExtra::combine(), Biobase::combine(), BiocGenerics::combine()
## x purrr::compact() masks XVector::compact()
## x purrr::compose() masks ShortRead::compose()
## x dplyr::count() masks matrixStats::count()
## x dplyr::desc() masks IRanges::desc()
## x tidyr::expand() masks S4Vectors::expand()
## x dplyr::filter() masks stats::filter()
## x dplyr::first() masks GenomicAlignments::first(), S4Vectors::first()
## x dplyr::id() masks ShortRead::id()
## x dplyr::lag() masks stats::lag()
## x dplyr::last() masks GenomicAlignments::last()
## x ggplot2::Position() masks BiocGenerics::Position(), base::Position()
## x purrr::reduce() masks GenomicRanges::reduce(), IRanges::reduce()
## x dplyr::rename() masks S4Vectors::rename()
## x purrr::simplify() masks DelayedArray::simplify()
## x dplyr::slice() masks XVector::slice(), IRanges::slice()
## x tibble::view() masks ShortRead::view()library(lubridate)##
## Attaching package: 'lubridate'## The following object is masked from 'package:GenomicAlignments':
##
## second## The following objects are masked from 'package:GenomicRanges':
##
## intersect, setdiff, union## The following object is masked from 'package:GenomeInfoDb':
##
## intersect## The following objects are masked from 'package:Biostrings':
##
## intersect, setdiff, union## The following objects are masked from 'package:IRanges':
##
## %within%, intersect, setdiff, union## The following objects are masked from 'package:S4Vectors':
##
## intersect, second, second<-, setdiff, union## The following objects are masked from 'package:BiocGenerics':
##
## intersect, setdiff, union## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(googlesheets)
library(here)## here() starts at /Users/marcel/OneDrive - Universität Zürich UZH/Experiment FK/multiple_stressors/multiple_stressorlibrary(scales)##
## Attaching package: 'scales'## The following object is masked from 'package:purrr':
##
## discard## The following object is masked from 'package:readr':
##
## col_factorlibrary(grid)
library(readxl)
library(ggplot2)
library(dplyr)
library(ggfortify)
library(ggpubr)
library(aod)
library(Rcpp)
library(dada2)
library(microbiome)##
## microbiome R package (microbiome.github.com)
##
##
##
## Copyright (C) 2011-2020 Leo Lahti,
## Sudarshan Shetty et al. <microbiome.github.io>##
## Attaching package: 'microbiome'## The following object is masked from 'package:scales':
##
## alpha## The following object is masked from 'package:ggplot2':
##
## alpha## The following object is masked from 'package:ShortRead':
##
## coverage## The following object is masked from 'package:GenomicAlignments':
##
## coverage## The following object is masked from 'package:SummarizedExperiment':
##
## coverage## The following object is masked from 'package:GenomicRanges':
##
## coverage## The following object is masked from 'package:Biostrings':
##
## coverage## The following objects are masked from 'package:IRanges':
##
## coverage, transform## The following object is masked from 'package:S4Vectors':
##
## transform## The following object is masked from 'package:base':
##
## transformlibrary(vegan)## Loading required package: permute## Loading required package: lattice##
## Attaching package: 'lattice'## The following object is masked from 'package:boot':
##
## melanoma## This is vegan 2.5-7##
## Attaching package: 'vegan'## The following object is masked from 'package:microbiome':
##
## diversitylibrary(plyr)## ------------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## ------------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following object is masked from 'package:ggpubr':
##
## mutate## The following object is masked from 'package:here':
##
## here## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize## The following object is masked from 'package:purrr':
##
## compact## The following object is masked from 'package:ShortRead':
##
## id## The following object is masked from 'package:matrixStats':
##
## count## The following object is masked from 'package:XVector':
##
## compact## The following object is masked from 'package:IRanges':
##
## desc## The following object is masked from 'package:S4Vectors':
##
## renamelibrary(DESeq2)
library(apeglm)
library(ggpubr)
library(patchwork)
library(broom)
library(rstatix)##
## Attaching package: 'rstatix'## The following objects are masked from 'package:plyr':
##
## desc, mutate## The following object is masked from 'package:IRanges':
##
## desc## The following object is masked from 'package:stats':
##
## filterlibrary(plyr)
library(permute)
library(lattice)
library(tm)## Loading required package: NLP##
## Attaching package: 'NLP'## The following object is masked from 'package:microbiome':
##
## meta## The following object is masked from 'package:ggplot2':
##
## annotate## The following objects are masked from 'package:Biobase':
##
## annotation, content## The following object is masked from 'package:BiocGenerics':
##
## annotation##
## Attaching package: 'tm'## The following object is masked from 'package:ShortRead':
##
## readerlibrary(stringr)
library(viridis) ## Loading required package: viridisLite##
## Attaching package: 'viridis'## The following object is masked from 'package:scales':
##
## viridis_palLoad phyloseq data, remove chloroplast sequences, remove samples that did not work
ps <- readRDS("phyloseq_multiple_stressor.rds")
taxa_names(ps) <- paste0("Seq", seq(ntaxa(ps)))
ps_new <- subset_taxa(ps, Order != "Chloroplast")
ps_new_1 <- subset_samples(ps_new,Lable != "ori_1")
ps_new_2 <- subset_samples(ps_new_1,Lable != "ori_4")
ps_new_3 <- subset_samples(ps_new_2,Lable != "B_1_24")
ps_new_4 <- subset_samples(ps_new_3,Lable != "C_1_20")
ps_new_5 <- subset_samples(ps_new_4,Lable != "BCD_5_24")
ps_new_6 <- subset_samples(ps_new_5,Lable != "CD_1_24")
ps_new_7 <- subset_samples(ps_new_6,Lable != "BE_5_24")
ps_new_8 <- subset_samples(ps_new_7,Lable != "BCDE_5_24")
ps_new_9 <- subset_samples(ps_new_8,Lable != "CD_5_24")
ps_new_10 <- subset_samples(ps_new_9,Lable != "BE_5_20")
ps_new_11 <- subset_samples(ps_new_10,Lable != "CE_1_24")Calculation of alpha diversity index
diversity_test<- estimate_richness(ps_new_11, split=TRUE, measures=NULL)## Warning in estimate_richness(ps_new_11, split = TRUE, measures = NULL): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.diversity_test$File <- rownames(diversity_test)
diversity_test$File <- gsub("\\.\\.", "--", diversity_test$File)
div_raw <- left_join(diversity_test,sample_data(ps_new_11)) ## Joining, by = "File"## Warning in class(x) <- c(setdiff(subclass, tibble_class), tibble_class): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer be an
## S4 objectdiv_raw <- div_raw %>%
reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))%>%
mutate(Temperature_factor=Temperature_factor)
#filter samples with low amount of reads
div_raw <- div_raw %>%
filter(Observed >=100) %>%
filter(Combination != "ori") %>%
mutate(T2=Temperature^2,T=Temperature_factor) %>%
select(-Temperature)%>%
reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))
stderr <- function(x, na.rm=FALSE) {
if (na.rm) x <- na.omit(x)
sqrt(var(x)/length(x))
}
#calculating mean of shannon and observed reads
div_means <- div_raw %>%
filter(T %in% c("20", "24", "28")) %>%
dplyr::group_by(Treatment, Combinations, T, Stressor,Stressor_numeric) %>%
dplyr::summarize(Shannon_mean = mean(Shannon),
Observed_mean = mean(Observed),
Shannon_SE = stderr(Shannon),
Observed_SE = stderr(Observed))%>%
reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))## `summarise()` has grouped output by 'Treatment', 'Combinations', 'T', 'Stressor'. You can override using the `.groups` argument.# plots
plot_diversity_all <- div_raw %>%
filter(T!="19")%>%
ggplot() +
geom_point(aes(x=Combinations, y=Shannon, shape=Stressor,col=Combinations), size = 3,position = position_dodge(width = 0.4)) +
facet_wrap("T")+
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold")) +
ylab ("Shannon index") +
xlab("Stressor(s) applied") +
scale_colour_viridis_d(direction = -1)+
theme(axis.text.x = element_text(angle = 50, hjust = 1))
plot_diversity_means <- div_means %>%
ggplot() +
geom_point(aes(x=T, y=Shannon_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=T,y=Shannon_mean,col=Combinations,ymin=Shannon_mean-Shannon_SE, ymax=Shannon_mean+Shannon_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
geom_line(aes(group=Combinations, x=T, y=Shannon_mean, col=Combinations),position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
scale_colour_viridis_d(direction = -1)+
xlab("Temperature [°C]") +
theme(legend.position = "none") +
ylab("Diversity_richness (mean)")
plot_diversity_means_numeric <- div_means %>%
ggplot() +
geom_point(aes(x=Stressor_numeric, y=Shannon_mean, col=T, shape= T), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=Shannon_mean,col=T,ymin=Shannon_mean-Shannon_SE, ymax=Shannon_mean+Shannon_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold")) +
xlab("Number of stressors") +
ylab("Diversity_richness (mean)") +
theme(legend.position = "none")
#final plot
(plot_diversity_means + plot_diversity_means_numeric)
#shannon on combinations
anova_diversity <- lm(Shannon ~F*G*M*A*T, data = div_raw)
autoplot(anova_diversity)## Warning: `arrange_()` was deprecated in dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
hist(residuals(anova_diversity ))
anova(anova_diversity )## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## F 1 4.155 4.1551 16.0171 9.316e-05 ***
## G 1 0.109 0.1093 0.4214 0.51709
## M 1 0.009 0.0090 0.0347 0.85238
## A 1 6.820 6.8197 26.2889 7.853e-07 ***
## T 2 3.539 1.7694 6.8208 0.00141 **
## F:G 1 0.200 0.2003 0.7721 0.38080
## F:M 1 0.055 0.0551 0.2125 0.64538
## G:M 1 0.076 0.0761 0.2935 0.58869
## F:A 1 4.638 4.6376 17.8770 3.820e-05 ***
## G:A 1 0.278 0.2778 1.0710 0.30216
## M:A 1 0.041 0.0408 0.1572 0.69228
## F:T 2 0.571 0.2853 1.0997 0.33531
## G:T 2 0.137 0.0684 0.2636 0.76857
## M:T 2 0.159 0.0797 0.3072 0.73594
## A:T 2 0.041 0.0204 0.0785 0.92455
## F:G:M 1 0.131 0.1313 0.5060 0.47784
## F:G:A 1 1.239 1.2386 4.7744 0.03024 *
## F:M:A 1 0.450 0.4500 1.7346 0.18958
## G:M:A 1 0.361 0.3614 1.3930 0.23953
## F:G:T 2 0.303 0.1516 0.5844 0.55852
## F:M:T 2 0.009 0.0044 0.0170 0.98311
## G:M:T 2 0.110 0.0548 0.2111 0.80990
## F:A:T 2 1.498 0.7491 2.8877 0.05841 .
## G:A:T 2 0.123 0.0616 0.2375 0.78884
## M:A:T 2 0.336 0.1678 0.6467 0.52502
## F:G:M:A 1 0.009 0.0087 0.0336 0.85481
## F:G:M:T 2 0.470 0.2349 0.9055 0.40628
## F:G:A:T 2 0.077 0.0387 0.1491 0.86158
## F:M:A:T 2 0.069 0.0347 0.1336 0.87502
## G:M:A:T 2 0.029 0.0143 0.0550 0.94646
## F:G:M:A:T 2 0.470 0.2352 0.9067 0.40577
## Residuals 172 44.619 0.2594
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(anova_diversity)##
## Call:
## lm(formula = Shannon ~ F * G * M * A * T, data = div_raw)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.93534 -0.20822 0.05585 0.27054 1.42868
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.521269 0.227778 24.240 <2e-16 ***
## F 0.074100 0.322127 0.230 0.8183
## G 0.364889 0.341668 1.068 0.2870
## M -0.040310 0.341668 -0.118 0.9062
## A 0.747909 0.341668 2.189 0.0299 *
## T24 0.146260 0.322127 0.454 0.6504
## T28 -0.180864 0.322127 -0.561 0.5752
## F:G -1.044270 0.469577 -2.224 0.0275 *
## F:M -0.109273 0.469577 -0.233 0.8163
## G:M -0.062624 0.496432 -0.126 0.8998
## F:A -0.526997 0.483191 -1.091 0.2769
## G:A -0.587326 0.483191 -1.216 0.2258
## M:A -0.166522 0.483191 -0.345 0.7308
## F:T24 -0.186145 0.469577 -0.396 0.6923
## F:T28 0.050948 0.455557 0.112 0.9111
## G:T24 -0.478917 0.483191 -0.991 0.3230
## G:T28 -0.393542 0.469577 -0.838 0.4032
## M:T24 -0.444855 0.469577 -0.947 0.3448
## M:T28 -0.263084 0.469577 -0.560 0.5760
## A:T24 -0.103835 0.483191 -0.215 0.8301
## A:T28 0.034482 0.469577 0.073 0.9415
## F:G:M 0.656103 0.673778 0.974 0.3315
## F:G:A 1.196845 0.683335 1.751 0.0816 .
## F:M:A 0.006661 0.673778 0.010 0.9921
## G:M:A 0.176000 0.683335 0.258 0.7971
## F:G:T24 1.390935 0.673778 2.064 0.0405 *
## F:G:T28 0.795634 0.654244 1.216 0.2256
## F:M:T24 0.852003 0.673778 1.265 0.2078
## F:M:T28 0.509086 0.654244 0.778 0.4376
## G:M:T24 0.837978 0.708192 1.183 0.2383
## G:M:T28 0.348712 0.673778 0.518 0.6054
## F:A:T24 -0.001744 0.692761 -0.003 0.9980
## F:A:T28 -0.493711 0.664083 -0.743 0.4582
## G:A:T24 0.571687 0.708192 0.807 0.4206
## G:A:T28 0.268085 0.664083 0.404 0.6869
## M:A:T24 0.759147 0.683335 1.111 0.2681
## M:A:T28 0.528232 0.664083 0.795 0.4275
## F:G:M:A -0.787214 0.952866 -0.826 0.4099
## F:G:M:T24 -1.818117 0.988499 -1.839 0.0676 .
## F:G:M:T28 -0.728523 0.932223 -0.781 0.4356
## F:G:A:T24 -1.259151 0.984115 -1.279 0.2025
## F:G:A:T28 -0.454171 0.939155 -0.484 0.6293
## F:M:A:T24 -0.650923 0.966382 -0.674 0.5015
## F:M:A:T28 -0.361584 0.932223 -0.388 0.6986
## G:M:A:T24 -1.157196 1.001535 -1.155 0.2495
## G:M:A:T28 -0.282198 0.939155 -0.300 0.7642
## F:G:M:A:T24 1.807982 1.394852 1.296 0.1966
## F:G:M:A:T28 0.425072 1.318363 0.322 0.7475
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5093 on 172 degrees of freedom
## Multiple R-squared: 0.3727, Adjusted R-squared: 0.2013
## F-statistic: 2.174 on 47 and 172 DF, p-value: 0.0001588Calculation rel abundance, just sequences that appear > 0.01 %, making dataframes for each column position
ps_rel <- transform_sample_counts(ps_new_11, function(x) x / sum(x) )
relab_threshold <- 0.001
ps_relab <- filter_taxa(ps_rel, function(x) !(sum(x < relab_threshold) == length(x)), TRUE)
ntaxa(ps_new_11)## [1] 20319ntaxa(ps_relab)## [1] 5830ps_relative <- transform_sample_counts(ps_relab, function(x) x / sum(x))NMDS
mds_whole <- ps_relative@otu_table %>%
as.data.frame() %>%
metaMDS(.,
distance = "bray",
k = 3, ## number of dimensions to reduce to
try = 100, ## number of random starts to try
autotransform = FALSE ## best not to use
)## Run 0 stress 0.1147651
## Run 1 stress 0.117328
## Run 2 stress 0.1201922
## Run 3 stress 0.1180983
## Run 4 stress 0.1171836
## Run 5 stress 0.11507
## ... Procrustes: rmse 0.006711343 max resid 0.08903122
## Run 6 stress 0.1192138
## Run 7 stress 0.1166306
## Run 8 stress 0.1223326
## Run 9 stress 0.1203552
## Run 10 stress 0.1148027
## ... Procrustes: rmse 0.002132423 max resid 0.01589763
## Run 11 stress 0.1158032
## Run 12 stress 0.1211516
## Run 13 stress 0.1147548
## ... New best solution
## ... Procrustes: rmse 0.002275816 max resid 0.02983052
## Run 14 stress 0.1165917
## Run 15 stress 0.1157718
## Run 16 stress 0.1221672
## Run 17 stress 0.1150052
## ... Procrustes: rmse 0.005425105 max resid 0.06951436
## Run 18 stress 0.1219735
## Run 19 stress 0.1149344
## ... Procrustes: rmse 0.00521367 max resid 0.0589545
## Run 20 stress 0.1157971
## Run 21 stress 0.1147687
## ... Procrustes: rmse 0.002162335 max resid 0.02918938
## Run 22 stress 0.1147639
## ... Procrustes: rmse 0.002125425 max resid 0.02817898
## Run 23 stress 0.1217448
## Run 24 stress 0.1162657
## Run 25 stress 0.1150722
## ... Procrustes: rmse 0.007000091 max resid 0.08936665
## Run 26 stress 0.1152504
## ... Procrustes: rmse 0.006393371 max resid 0.07191682
## Run 27 stress 0.1147718
## ... Procrustes: rmse 0.002082273 max resid 0.0268729
## Run 28 stress 0.1147541
## ... New best solution
## ... Procrustes: rmse 0.0002902666 max resid 0.003265375
## ... Similar to previous best
## Run 29 stress 0.1151349
## ... Procrustes: rmse 0.006531455 max resid 0.06005487
## Run 30 stress 0.1166066
## Run 31 stress 0.1147726
## ... Procrustes: rmse 0.001761228 max resid 0.02059767
## Run 32 stress 0.1157936
## Run 33 stress 0.1148025
## ... Procrustes: rmse 0.004203744 max resid 0.05364829
## Run 34 stress 0.1166647
## Run 35 stress 0.1178293
## Run 36 stress 0.1151829
## ... Procrustes: rmse 0.006275269 max resid 0.06002028
## Run 37 stress 0.1150307
## ... Procrustes: rmse 0.00493667 max resid 0.04736456
## Run 38 stress 0.1165166
## Run 39 stress 0.1149577
## ... Procrustes: rmse 0.00367769 max resid 0.03573953
## Run 40 stress 0.1186977
## Run 41 stress 0.1237333
## Run 42 stress 0.1164721
## Run 43 stress 0.1150399
## ... Procrustes: rmse 0.006451346 max resid 0.08959954
## Run 44 stress 0.1272646
## Run 45 stress 0.1271736
## Run 46 stress 0.1148437
## ... Procrustes: rmse 0.003686582 max resid 0.04045134
## Run 47 stress 0.1149457
## ... Procrustes: rmse 0.004186228 max resid 0.04788836
## Run 48 stress 0.1215899
## Run 49 stress 0.1216396
## Run 50 stress 0.1161779
## Run 51 stress 0.1147746
## ... Procrustes: rmse 0.002987398 max resid 0.04055853
## Run 52 stress 0.1233454
## Run 53 stress 0.1232819
## Run 54 stress 0.1160627
## Run 55 stress 0.1148072
## ... Procrustes: rmse 0.002016763 max resid 0.02421261
## Run 56 stress 0.1147672
## ... Procrustes: rmse 0.001491363 max resid 0.01096601
## Run 57 stress 0.1149019
## ... Procrustes: rmse 0.004133279 max resid 0.03957484
## Run 58 stress 0.1148842
## ... Procrustes: rmse 0.003865038 max resid 0.03467656
## Run 59 stress 0.1147631
## ... Procrustes: rmse 0.00103653 max resid 0.007984618
## ... Similar to previous best
## Run 60 stress 0.1158119
## Run 61 stress 0.120041
## Run 62 stress 0.1216382
## Run 63 stress 0.1151596
## ... Procrustes: rmse 0.006141564 max resid 0.08723947
## Run 64 stress 0.1149327
## ... Procrustes: rmse 0.00426985 max resid 0.05363242
## Run 65 stress 0.12263
## Run 66 stress 0.1227652
## Run 67 stress 0.1219891
## Run 68 stress 0.1163769
## Run 69 stress 0.114763
## ... Procrustes: rmse 0.00106261 max resid 0.0112492
## Run 70 stress 0.1149026
## ... Procrustes: rmse 0.005905025 max resid 0.05642366
## Run 71 stress 0.1190239
## Run 72 stress 0.116873
## Run 73 stress 0.11577
## Run 74 stress 0.1206862
## Run 75 stress 0.1196646
## Run 76 stress 0.1147574
## ... Procrustes: rmse 0.0006005219 max resid 0.004559622
## ... Similar to previous best
## Run 77 stress 0.1168291
## Run 78 stress 0.1224913
## Run 79 stress 0.1161734
## Run 80 stress 0.1227812
## Run 81 stress 0.1164785
## Run 82 stress 0.1241439
## Run 83 stress 0.1165355
## Run 84 stress 0.1155676
## Run 85 stress 0.1158396
## Run 86 stress 0.1213225
## Run 87 stress 0.1148268
## ... Procrustes: rmse 0.002679179 max resid 0.02778069
## Run 88 stress 0.1148036
## ... Procrustes: rmse 0.002636979 max resid 0.03015428
## Run 89 stress 0.1196198
## Run 90 stress 0.1157682
## Run 91 stress 0.1254822
## Run 92 stress 0.1179103
## Run 93 stress 0.1147765
## ... Procrustes: rmse 0.001213848 max resid 0.00936875
## ... Similar to previous best
## Run 94 stress 0.1148853
## ... Procrustes: rmse 0.004558338 max resid 0.0541002
## Run 95 stress 0.1169937
## Run 96 stress 0.1230001
## Run 97 stress 0.1201933
## Run 98 stress 0.1168677
## Run 99 stress 0.1221757
## Run 100 stress 0.1271232
## *** Solution reached ## 0.11
mds_whole_res <- ps_relative@sam_data %>%
as.tibble() %>%
select(Treatment, Stressor,Stressor_numeric, Temperature, Combinations, Temperature_factor, Lable, Name,F,G,M,A) %>%
bind_cols(as.tibble(scores(mds_whole, display = "sites"))) %>%
mutate(T2=Temperature^2,T=Temperature_factor) %>%
select(-Temperature)## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.## Warning in class(x) <- c(setdiff(subclass, tibble_class), tibble_class): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer be an
## S4 objectmds_whole_res <- mds_whole_res %>%
filter(Temperature_factor != "19") %>%
reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))
ggplot(mds_whole_res, aes(x = NMDS1, y = NMDS2)) +
geom_point(size = 2) +
facet_grid(Combinations~Temperature_factor)+
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))
Plots of NMDS distances that show high significance
#means of NMDS
means_NMDS <- mds_whole_res %>%
filter(Temperature_factor %in% c("20", "24", "28")) %>%
dplyr::group_by(Treatment,Combinations, Stressor, Stressor_numeric, Temperature_factor,F,G,M,A) %>%
dplyr::summarize(NMDS1_mean = mean(NMDS1),
NMDS2_mean = mean(NMDS2),
NMDS3_mean = mean(NMDS3),
NMDS1_SE = stderr(NMDS1),
NMDS2_SE = stderr(NMDS2),
NMDS3_SE= stderr(NMDS3))## `summarise()` has grouped output by 'Treatment', 'Combinations', 'Stressor', 'Stressor_numeric', 'Temperature_factor', 'F', 'G', 'M'. You can override using the `.groups` argument.NMDS1 and NMDS2 plots
#NMDS1 supplement plot
nmds1_plot_temperature <- means_NMDS %>%
ggplot() +
geom_point(aes(x=Temperature_factor, y=NMDS1_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Temperature_factor,y=NMDS1_mean,col=Combinations,ymin=NMDS1_mean-NMDS1_SE, ymax=NMDS1_mean+NMDS1_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
geom_line(aes(group=Combinations, x=Temperature_factor, y=NMDS1_mean, col=Combinations),position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
scale_colour_viridis_d(direction = -1)+
xlab("Temperature") +
ylab("NMDS1 (mean)") +
theme(legend.position = "none")
plot_NMDS1_means_numeric <- means_NMDS %>%
ggplot() +
geom_point(aes(x=Stressor_numeric, y=NMDS1_mean, col=Temperature_factor, shape= Temperature_factor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=NMDS1_mean,col=Temperature_factor,ymin=NMDS1_mean-NMDS1_SE, ymax=NMDS1_mean+NMDS1_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
xlab("Number of stressors") +
ylab("NMDS1 (mean)") +
theme(legend.position = "none")
# plot NMDS1
nmds1_plot_temperature + plot_NMDS1_means_numeric 
#NMDS2 supplement plot
NMDS2_plot_temperature <- means_NMDS %>%
ggplot() +
geom_point(aes(x=Temperature_factor, y=NMDS2_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Temperature_factor,y=NMDS2_mean,col=Combinations,ymin=NMDS2_mean-NMDS2_SE, ymax=NMDS2_mean+NMDS2_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
geom_line(aes(group=Combinations, x=Temperature_factor, y=NMDS2_mean, col=Combinations),position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
scale_colour_viridis_d(direction = -1)+
xlab("Temperature") +
theme(legend.position = "none") +
ylab("NMDS2 (mean)")
plot_NMDS2_means_numeric <- means_NMDS %>%
ggplot() +
geom_point(aes(x=Stressor_numeric, y=NMDS2_mean, col=Temperature_factor, shape= Temperature_factor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=NMDS2_mean,col=Temperature_factor,ymin=NMDS2_mean-NMDS2_SE, ymax=NMDS2_mean+NMDS2_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold")) +
xlab("Number of stressors") +
ylab("NMDS2 (mean)")+
theme(legend.position = "none")
# plot NMDS2
NMDS2_plot_temperature + plot_NMDS2_means_numeric 
###NMDS1 anova and plots
anova_NMDS1 <- lm(NMDS1 ~F*G*M*A*T, data = mds_whole_res)
autoplot(anova_NMDS1)
hist(residuals(anova_NMDS1 ))
anova(anova_NMDS1 )## Analysis of Variance Table
##
## Response: NMDS1
## Df Sum Sq Mean Sq F value Pr(>F)
## F 1 61.036 61.036 1028.3184 < 2.2e-16 ***
## G 1 0.002 0.002 0.0393 0.843076
## M 1 0.029 0.029 0.4889 0.485297
## A 1 2.625 2.625 44.2292 3.266e-10 ***
## T 2 0.652 0.326 5.4883 0.004844 **
## F:G 1 0.220 0.220 3.7031 0.055862 .
## F:M 1 0.000 0.000 0.0009 0.976463
## G:M 1 0.142 0.142 2.4002 0.123044
## F:A 1 3.009 3.009 50.7008 2.375e-11 ***
## G:A 1 0.014 0.014 0.2341 0.629091
## M:A 1 0.023 0.023 0.3956 0.530179
## F:T 2 2.659 1.329 22.3970 1.991e-09 ***
## G:T 2 0.068 0.034 0.5768 0.562723
## M:T 2 0.063 0.031 0.5270 0.591274
## A:T 2 0.211 0.105 1.7773 0.172005
## F:G:M 1 0.001 0.001 0.0203 0.886743
## F:G:A 1 0.128 0.128 2.1539 0.143927
## F:M:A 1 0.011 0.011 0.1884 0.664797
## G:M:A 1 0.130 0.130 2.1874 0.140864
## F:G:T 2 0.100 0.050 0.8402 0.433284
## F:M:T 2 0.033 0.017 0.2821 0.754527
## G:M:T 2 0.423 0.212 3.5659 0.030251 *
## F:A:T 2 0.452 0.226 3.8116 0.023888 *
## G:A:T 2 0.125 0.063 1.0551 0.350268
## M:A:T 2 0.037 0.019 0.3119 0.732470
## F:G:M:A 1 0.001 0.001 0.0243 0.876249
## F:G:M:T 2 0.046 0.023 0.3859 0.680386
## F:G:A:T 2 0.030 0.015 0.2526 0.777054
## F:M:A:T 2 0.016 0.008 0.1310 0.877296
## G:M:A:T 2 0.181 0.090 1.5243 0.220517
## F:G:M:A:T 2 0.021 0.011 0.1779 0.837146
## Residuals 183 10.862 0.059
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(anova_NMDS1)##
## Call:
## lm(formula = NMDS1 ~ F * G * M * A * T, data = mds_whole_res)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.90640 -0.12562 -0.00319 0.08877 0.98663
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.62357 0.10895 -5.723 4.21e-08 ***
## F 1.01436 0.15408 6.583 4.72e-10 ***
## G 0.08727 0.16343 0.534 0.5940
## M -0.10535 0.15408 -0.684 0.4950
## A 0.25032 0.15408 1.625 0.1060
## T24 0.04710 0.15408 0.306 0.7602
## T28 -0.15322 0.15408 -0.994 0.3213
## F:G -0.03967 0.22462 -0.177 0.8600
## F:M 0.15503 0.21791 0.711 0.4777
## G:M -0.07758 0.22462 -0.345 0.7302
## F:A -0.34377 0.22462 -1.531 0.1276
## G:A -0.08645 0.22462 -0.385 0.7008
## M:A 0.07985 0.21791 0.366 0.7145
## F:T24 0.11684 0.22462 0.520 0.6036
## F:T28 0.51957 0.21791 2.384 0.0181 *
## G:T24 -0.40315 0.22462 -1.795 0.0743 .
## G:T28 -0.08559 0.22462 -0.381 0.7036
## M:T24 -0.11534 0.21791 -0.529 0.5972
## M:T28 -0.02001 0.21791 -0.092 0.9269
## A:T24 0.30010 0.21791 1.377 0.1701
## A:T28 0.20718 0.21791 0.951 0.3430
## F:G:M -0.03538 0.31295 -0.113 0.9101
## F:G:A 0.32514 0.31765 1.024 0.3074
## F:M:A -0.12032 0.31295 -0.384 0.7011
## G:M:A 0.01457 0.31295 0.047 0.9629
## F:G:T24 0.19467 0.31765 0.613 0.5407
## F:G:T28 -0.01761 0.31295 -0.056 0.9552
## F:M:T24 -0.26511 0.31295 -0.847 0.3980
## F:M:T28 -0.11238 0.30817 -0.365 0.7158
## G:M:T24 0.62782 0.32535 1.930 0.0552 .
## G:M:T28 0.16828 0.31295 0.538 0.5914
## F:A:T24 -0.36116 0.32229 -1.121 0.2639
## F:A:T28 -0.20663 0.31295 -0.660 0.5099
## G:A:T24 0.18880 0.31765 0.594 0.5530
## G:A:T28 -0.05402 0.31295 -0.173 0.8631
## M:A:T24 0.20496 0.30817 0.665 0.5068
## M:A:T28 0.02853 0.30817 0.093 0.9263
## F:G:M:A -0.07491 0.44257 -0.169 0.8658
## F:G:M:T24 0.09541 0.45470 0.210 0.8340
## F:G:M:T28 0.03307 0.43921 0.075 0.9401
## F:G:A:T24 -0.33135 0.45252 -0.732 0.4650
## F:G:A:T28 -0.01474 0.44257 -0.033 0.9735
## F:M:A:T24 0.06044 0.44592 0.136 0.8923
## F:M:A:T28 0.19036 0.43921 0.433 0.6652
## G:M:A:T24 -0.62681 0.45143 -1.388 0.1667
## G:M:A:T28 0.02288 0.43921 0.052 0.9585
## F:G:M:A:T24 0.24337 0.63919 0.381 0.7038
## F:G:M:A:T28 -0.13143 0.62114 -0.212 0.8327
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2436 on 183 degrees of freedom
## Multiple R-squared: 0.8697, Adjusted R-squared: 0.8362
## F-statistic: 25.98 on 47 and 183 DF, p-value: < 2.2e-16#F:A:Temperature_factor
ggplot(means_NMDS) +
geom_line(aes(x=Temperature_factor,y=NMDS1_mean, col=as.factor(F),shape=as.factor(A),
group=interaction(as.factor(F),A,M,G)), alpha=0.2)+
geom_point(mapping = aes(x=Temperature_factor,y=NMDS1_mean, col=as.factor(F),shape=as.factor(A)),size=3)+
theme_bw() + theme(panel.grid = element_blank())+
geom_errorbar(aes(x=Temperature_factor,y=NMDS1_mean, ymin=NMDS1_mean-NMDS1_SE,ymax=NMDS1_mean+NMDS1_SE,col=as.factor(F)),width=0.1,alpha=0.2)+
xlab("Temperature [°C]")## Warning: Ignoring unknown aesthetics: shape
###NMDS2 anova and plots
#NMDS2 combination
anova_NMDS2 <- lm(NMDS2 ~F*G*M*A*T, data = mds_whole_res)
autoplot(anova_NMDS2)
hist(residuals(anova_NMDS2 ))
anova(anova_NMDS2 )## Analysis of Variance Table
##
## Response: NMDS2
## Df Sum Sq Mean Sq F value Pr(>F)
## F 1 0.0493 0.0493 0.3960 0.529963
## G 1 0.0004 0.0004 0.0034 0.953514
## M 1 0.7552 0.7552 6.0673 0.014695 *
## A 1 4.8989 4.8989 39.3559 2.484e-09 ***
## T 2 7.7550 3.8775 31.1500 2.275e-12 ***
## F:G 1 0.1971 0.1971 1.5831 0.209918
## F:M 1 0.0523 0.0523 0.4200 0.517749
## G:M 1 0.0059 0.0059 0.0471 0.828501
## F:A 1 0.9378 0.9378 7.5338 0.006659 **
## G:A 1 0.4115 0.4115 3.3061 0.070658 .
## M:A 1 0.1714 0.1714 1.3766 0.242201
## F:T 2 0.0725 0.0362 0.2911 0.747764
## G:T 2 0.2108 0.1054 0.8467 0.430497
## M:T 2 0.2787 0.1394 1.1196 0.328627
## A:T 2 0.2865 0.1433 1.1508 0.318652
## F:G:M 1 0.0888 0.0888 0.7132 0.399501
## F:G:A 1 0.0137 0.0137 0.1101 0.740426
## F:M:A 1 0.3646 0.3646 2.9292 0.088686 .
## G:M:A 1 0.4031 0.4031 3.2382 0.073587 .
## F:G:T 2 0.0916 0.0458 0.3678 0.692757
## F:M:T 2 0.1976 0.0988 0.7937 0.453731
## G:M:T 2 0.3708 0.1854 1.4896 0.228179
## F:A:T 2 1.1742 0.5871 4.7163 0.010064 *
## G:A:T 2 0.4435 0.2217 1.7813 0.171333
## M:A:T 2 0.1121 0.0561 0.4503 0.638123
## F:G:M:A 1 0.0078 0.0078 0.0624 0.803066
## F:G:M:T 2 0.0597 0.0299 0.2399 0.786919
## F:G:A:T 2 0.2561 0.1281 1.0287 0.359528
## F:M:A:T 2 0.0270 0.0135 0.1084 0.897316
## G:M:A:T 2 0.1251 0.0625 0.5025 0.605855
## F:G:M:A:T 2 0.7768 0.3884 3.1202 0.046506 *
## Residuals 183 22.7795 0.1245
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(anova_NMDS2)##
## Call:
## lm(formula = NMDS2 ~ F * G * M * A * T, data = mds_whole_res)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.67576 -0.10981 -0.00263 0.11760 1.33865
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.159311 0.157783 -1.010 0.3140
## F 0.085881 0.223139 0.385 0.7008
## G 0.067201 0.236675 0.284 0.7768
## M 0.126427 0.223139 0.567 0.5717
## A -0.231951 0.223139 -1.039 0.2999
## T24 0.434002 0.223139 1.945 0.0533 .
## T28 0.554193 0.223139 2.484 0.0139 *
## F:G 0.393576 0.325279 1.210 0.2279
## F:M -0.274407 0.315567 -0.870 0.3857
## G:M -0.153747 0.325279 -0.473 0.6370
## F:A -0.038816 0.325279 -0.119 0.9051
## G:A 0.076456 0.325279 0.235 0.8144
## M:A -0.027077 0.315567 -0.086 0.9317
## F:T24 -0.103381 0.325279 -0.318 0.7510
## F:T28 -0.230081 0.315567 -0.729 0.4669
## G:T24 0.144518 0.325279 0.444 0.6574
## G:T28 -0.032296 0.325279 -0.099 0.9210
## M:T24 -0.144553 0.315567 -0.458 0.6474
## M:T28 -0.088370 0.315567 -0.280 0.7798
## A:T24 0.094618 0.315567 0.300 0.7646
## A:T28 -0.167913 0.315567 -0.532 0.5953
## F:G:M -0.335092 0.453198 -0.739 0.4606
## F:G:A -0.647511 0.460014 -1.408 0.1609
## F:M:A 0.149664 0.453198 0.330 0.7416
## G:M:A -0.017552 0.453198 -0.039 0.9691
## F:G:T24 -0.742368 0.460014 -1.614 0.1083
## F:G:T28 -0.181871 0.453198 -0.401 0.6887
## F:M:T24 -0.227877 0.453198 -0.503 0.6157
## F:M:T28 0.192243 0.446279 0.431 0.6671
## G:M:T24 -0.097784 0.471154 -0.208 0.8358
## G:M:T28 0.031290 0.453198 0.069 0.9450
## F:A:T24 -0.009746 0.466730 -0.021 0.9834
## F:A:T28 0.541322 0.453198 1.194 0.2338
## G:A:T24 -1.011734 0.460014 -2.199 0.0291 *
## G:A:T28 -0.090973 0.453198 -0.201 0.8411
## M:A:T24 -0.449489 0.446279 -1.007 0.3152
## M:A:T28 -0.044948 0.446279 -0.101 0.9199
## F:G:M:A 0.580079 0.640919 0.905 0.3666
## F:G:M:T24 1.183700 0.658482 1.798 0.0739 .
## F:G:M:T28 0.151148 0.636045 0.238 0.8124
## F:G:A:T24 1.576498 0.655324 2.406 0.0171 *
## F:G:A:T28 0.138070 0.640919 0.215 0.8297
## F:M:A:T24 0.691480 0.645756 1.071 0.2857
## F:M:A:T28 -0.283331 0.636045 -0.445 0.6565
## G:M:A:T24 0.766263 0.653739 1.172 0.2427
## G:M:A:T28 0.153726 0.636045 0.242 0.8093
## F:G:M:A:T24 -1.853952 0.925648 -2.003 0.0467 *
## F:G:M:A:T28 0.279283 0.899504 0.310 0.7565
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3528 on 183 degrees of freedom
## Multiple R-squared: 0.4748, Adjusted R-squared: 0.3399
## F-statistic: 3.52 on 47 and 183 DF, p-value: 6.71e-10Visulaization of relative abundance of family level of all combinations with heatmap
df_rel <- psmelt(ps_relative)
df_family<- df_rel %>%
group_by(Family,Lable,Combinations,Temperature_factor) %>%
dplyr::summarize(Abundance=sum(Abundance))## `summarise()` has grouped output by 'Family', 'Lable', 'Combinations'. You can override using the `.groups` argument.df_family <- df_family %>%
reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))
df_family%>%
filter(Family %in% c("Phormidiaceae","Chromatiaceae"),Combinations != "ori") %>%
ggplot(aes(x=Combinations, y=Abundance, fill=Family)) +
geom_boxplot()+
facet_grid("Temperature_factor") +
theme_bw()
df_family <- df_family %>%
dplyr::group_by(Family,Combinations,Temperature_factor) %>%
dplyr::summarize(Abundance_mean = mean(Abundance),
Abundance_SE =stderr(Abundance))## `summarise()` has grouped output by 'Family', 'Combinations'. You can override using the `.groups` argument.index <- which(df_family$Abundance_mean>=0.05)
family_to_keep <- unique(df_family[index,"Family"])
family_to_keep <- unname(unlist(family_to_keep))
df_family $Family_filter <- ifelse(df_family$Family %in% family_to_keep, df_family$Family,"other")
#make other to one
df_family<- df_family %>%
dplyr::group_by(Family_filter,Combinations,Temperature_factor) %>%
dplyr::summarize(Abundance_sum = sum(Abundance_mean),
Abundance_sum_SE =sum(Abundance_SE))## `summarise()` has grouped output by 'Family_filter', 'Combinations'. You can override using the `.groups` argument.df_family <- df_family %>%
reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))
cols <-c("Azospirillaceae"="bisque1", "Bacteroidetes_vadinHA17" = "yellow", "Chlorobiaceae"="black","Comamonadaceae" = "green4", "Hungateiclostridiaceae" ="maroon1", "Lentimicrobiaceae" ="royalblue1","Marinifilaceae"="lightcyan1","Nocardiaceae" ="thistle1", "Phormidiaceae"= "firebrick1", "Prolixibacteraceae" = "magenta4", "Pseudomonadaceae" = "grey","Rhodocyclaceae"="green", "Rhodospillaceae"="grey38","Rubritaleaceae"="brown", "Xanthobacteraceae"="cornflowerblue", "Xanthomonadaceae"="lightskyblue1","Cyanobiaceae"="lavender", "Chromatiaceae" = "lightsalmon")
community_plot <- df_family %>%
filter(Combinations != "ori",Family_filter !="Mitochondria",Family_filter !="other")%>%
ggplot(aes(x = Combinations, y = Abundance_sum, fill=Family_filter)) +
geom_col() +
facet_wrap("Temperature_factor")+
theme_bw()+
theme(axis.text.x = element_text(angle = 50, hjust = 1,size=14)) +
scale_fill_manual("Family", values = cols) +
ylab("rel. abundance family level") +
xlab("Stressor(s) applied") +
theme(axis.title=element_text(size=16,face="bold")) +
theme(legend.text = element_text(size=11))
df_family %>%
filter(Combinations != "ori", Family_filter != "other", Family_filter !="Pirellulaceae", Family_filter !="Cyanobiaceae",Family_filter !="Mitochondria")%>%
ggplot(aes(x = Combinations, y = Family_filter, fill=Abundance_sum)) +
geom_tile() +
scale_fill_gradient2(low="red", mid="white", high="navy", midpoint=0.15, na.value="lightblue",
breaks=c(0.05,0.1,0.15,0.2,0.25),
limits=c(0.05,0.25)) +
facet_wrap("Temperature_factor")+
theme(strip.placement = "inside") +
theme_bw()+
theme(axis.text.x = element_text(angle = 50, hjust = 1))
#just phormidiaceae and chromatiaceae for supplementary information
df_family %>%
filter(Family_filter %in% c("Phormidiaceae","Chromatiaceae"), Combinations!= "ori") %>%
ggplot() +
geom_point(aes(x=Combinations,y=Abundance_sum,col=Family_filter,fill=Family_filter),position=position_dodge(width = 0.1),size=4)+
geom_errorbar(aes(x=Combinations,y=Abundance_sum,col=Family_filter,ymin=Abundance_sum-Abundance_sum_SE, ymax=Abundance_sum+Abundance_sum_SE),width=0.1,alpha=0.4,col="black", position=position_dodge(width = .1))+
facet_wrap("Temperature_factor")+
theme_bw() +
ylab("rel. abundance") +
theme(axis.text.x = element_text(angle = 50, hjust = 1, size=14)) +
labs(col = "Family") + labs(fill="Family")+
xlab("Stressor combination")
Genus visualization
df_genus<- df_rel %>%
group_by(Genus,Lable,Combinations,Temperature_factor) %>%
dplyr::summarize(Abundance=sum(Abundance))## `summarise()` has grouped output by 'Genus', 'Lable', 'Combinations'. You can override using the `.groups` argument.df_genus <- df_genus %>%
dplyr::group_by(Genus,Combinations,Temperature_factor) %>%
dplyr::summarize(Abundance = mean(Abundance))## `summarise()` has grouped output by 'Genus', 'Combinations'. You can override using the `.groups` argument.index <- which(df_genus$Abundance>=0.05)
genus_to_keep <- unique(df_genus[index,"Genus"])
genus_to_keep <- unname(unlist(genus_to_keep))
df_genus $Genus_filter <- ifelse(df_genus$Genus %in% genus_to_keep, df_genus$Genus,"other")
#make other to one
df_genus<- df_genus %>%
dplyr::group_by(Genus_filter,Combinations,Temperature_factor) %>%
dplyr::summarize(Abundance = sum(Abundance))## `summarise()` has grouped output by 'Genus_filter', 'Combinations'. You can override using the `.groups` argument.df_genus <- df_genus %>%
reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))
genus_plot <- df_genus %>%
filter(Combinations != "ori", Genus_filter != "other", Genus_filter !="Pirellulaceae", Genus_filter !="Cyanobium_PCC-6307", Genus_filter !="HN-HF0106")%>%
ggplot(aes(x = Combinations, y = Genus_filter, fill=Abundance)) +
geom_tile() +
scale_fill_continuous(name = "rel.abundance",
low = "#FF0000",
high = "#00CCFF",
na.value="white",
breaks=c(0.05,0.1,0.15,0.2,0.25),
limits=c(0.05,0.25)) +
facet_wrap("Temperature_factor") +
theme(axis.text.x = element_text(angle = 50, hjust = 1, size=14)) +
theme(axis.text.y = element_text(size=12,face="italic"))+
theme(axis.title=element_text(size=16,face="bold")) +
xlab("Stressor(s) applied") +
ylab("rel. abundance genus level")(community_plot) +
(plot_spacer() + genus_plot + plot_layout(widths=c(1,100))) +
plot_layout(ncol=1,heights=c(1,1)) + plot_annotation(tag_levels = "a")
Subpupulation dataset
ps_rel <- transform_sample_counts(ps_new_11, function(x) x / sum(x) )
relab_threshold <- 0.001
ps_relab <- filter_taxa(ps_rel, function(x) !(sum(x < relab_threshold) == length(x)), TRUE)
ntaxa(ps_new)## [1] 20319ntaxa(ps_relab)## [1] 5830ps_interest <- subset_taxa(ps_relab, Order %in% c("Cyanobacteriales","Chromatiales","Synechococcales", "Chlorobiales","Leptolyngbyales","Desulfobulbales","Desulfovibrionales","Limnotrichales","Campylobacterales"))
ntaxa(ps_interest)## [1] 493ps_relative_interest <- transform_sample_counts(ps_interest, function(x) x / sum(x))
interest <-psmelt(ps_relative_interest)
ps_nmds <- subset_samples(ps_relative_interest, Treatment != "ori")
mds_whole_subset <- ps_nmds@otu_table %>%
as.data.frame() %>%
metaMDS(.,
distance = "bray",
k = 3, ## number of dimensions to reduce to
try = 300, ## number of random starts to try
autotransform = FALSE ## best not to use
) ## Run 0 stress 0.1048781
## Run 1 stress 0.1072036
## Run 2 stress 0.1037414
## ... New best solution
## ... Procrustes: rmse 0.006635281 max resid 0.08474109
## Run 3 stress 0.105289
## Run 4 stress 0.1051615
## Run 5 stress 0.1061844
## Run 6 stress 0.1063321
## Run 7 stress 0.1072051
## Run 8 stress 0.1049566
## Run 9 stress 0.1052647
## Run 10 stress 0.106756
## Run 11 stress 0.1037896
## ... Procrustes: rmse 0.00293663 max resid 0.02417046
## Run 12 stress 0.1065792
## Run 13 stress 0.1071902
## Run 14 stress 0.1056136
## Run 15 stress 0.1073591
## Run 16 stress 0.1057231
## Run 17 stress 0.1052825
## Run 18 stress 0.1064098
## Run 19 stress 0.1053836
## Run 20 stress 0.1048949
## Run 21 stress 0.1048907
## Run 22 stress 0.1052388
## Run 23 stress 0.1051941
## Run 24 stress 0.1049397
## Run 25 stress 0.105093
## Run 26 stress 0.1037102
## ... New best solution
## ... Procrustes: rmse 0.003714632 max resid 0.02826341
## Run 27 stress 0.1054493
## Run 28 stress 0.1053349
## Run 29 stress 0.1037115
## ... Procrustes: rmse 0.0005342914 max resid 0.007379712
## ... Similar to previous best
## Run 30 stress 0.1039884
## ... Procrustes: rmse 0.008394872 max resid 0.1217423
## Run 31 stress 0.103687
## ... New best solution
## ... Procrustes: rmse 0.001862688 max resid 0.02555876
## Run 32 stress 0.1040644
## ... Procrustes: rmse 0.008425929 max resid 0.1189921
## Run 33 stress 0.1037651
## ... Procrustes: rmse 0.002969157 max resid 0.04276497
## Run 34 stress 0.1052573
## Run 35 stress 0.1066791
## Run 36 stress 0.1036805
## ... New best solution
## ... Procrustes: rmse 0.0007859608 max resid 0.009224061
## ... Similar to previous best
## Run 37 stress 0.1066239
## Run 38 stress 0.1051974
## Run 39 stress 0.1072624
## Run 40 stress 0.103879
## ... Procrustes: rmse 0.004549549 max resid 0.043179
## Run 41 stress 0.1039829
## ... Procrustes: rmse 0.007968812 max resid 0.117865
## Run 42 stress 0.1054307
## Run 43 stress 0.1065277
## Run 44 stress 0.1063728
## Run 45 stress 0.1071685
## Run 46 stress 0.1053852
## Run 47 stress 0.1051288
## Run 48 stress 0.1072725
## Run 49 stress 0.1055484
## Run 50 stress 0.1054812
## Run 51 stress 0.106358
## Run 52 stress 0.1053356
## Run 53 stress 0.1039857
## ... Procrustes: rmse 0.008051974 max resid 0.1201322
## Run 54 stress 0.1063092
## Run 55 stress 0.1037752
## ... Procrustes: rmse 0.002983071 max resid 0.04208786
## Run 56 stress 0.1040888
## ... Procrustes: rmse 0.008837648 max resid 0.1238711
## Run 57 stress 0.1051795
## Run 58 stress 0.1053434
## Run 59 stress 0.1038061
## ... Procrustes: rmse 0.003754823 max resid 0.0245486
## Run 60 stress 0.1040709
## ... Procrustes: rmse 0.007758431 max resid 0.1073494
## Run 61 stress 0.1052072
## Run 62 stress 0.107323
## Run 63 stress 0.1051288
## Run 64 stress 0.1037687
## ... Procrustes: rmse 0.002974885 max resid 0.04227947
## Run 65 stress 0.1037411
## ... Procrustes: rmse 0.003008803 max resid 0.02138887
## Run 66 stress 0.1067213
## Run 67 stress 0.1052758
## Run 68 stress 0.1052533
## Run 69 stress 0.1096792
## Run 70 stress 0.1054532
## Run 71 stress 0.1051643
## Run 72 stress 0.105183
## Run 73 stress 0.1040851
## ... Procrustes: rmse 0.008028575 max resid 0.1092318
## Run 74 stress 0.1037936
## ... Procrustes: rmse 0.003547532 max resid 0.04232409
## Run 75 stress 0.1049463
## Run 76 stress 0.1065428
## Run 77 stress 0.1037843
## ... Procrustes: rmse 0.003296965 max resid 0.04294894
## Run 78 stress 0.1054323
## Run 79 stress 0.1054745
## Run 80 stress 0.1039897
## ... Procrustes: rmse 0.007517483 max resid 0.1102222
## Run 81 stress 0.1038034
## ... Procrustes: rmse 0.003542826 max resid 0.02124059
## Run 82 stress 0.1037411
## ... Procrustes: rmse 0.00300768 max resid 0.02143947
## Run 83 stress 0.1051659
## Run 84 stress 0.1038814
## ... Procrustes: rmse 0.004661334 max resid 0.04320149
## Run 85 stress 0.1040868
## ... Procrustes: rmse 0.008775891 max resid 0.1229213
## Run 86 stress 0.1054308
## Run 87 stress 0.1039924
## ... Procrustes: rmse 0.008204912 max resid 0.122926
## Run 88 stress 0.1036944
## ... Procrustes: rmse 0.001268815 max resid 0.01526923
## Run 89 stress 0.1052637
## Run 90 stress 0.1055796
## Run 91 stress 0.1052105
## Run 92 stress 0.1040023
## ... Procrustes: rmse 0.004446841 max resid 0.04310146
## Run 93 stress 0.1040637
## ... Procrustes: rmse 0.00846594 max resid 0.1184098
## Run 94 stress 0.1039327
## ... Procrustes: rmse 0.005352253 max resid 0.06325116
## Run 95 stress 0.1044426
## Run 96 stress 0.1052806
## Run 97 stress 0.1054762
## Run 98 stress 0.1065128
## Run 99 stress 0.1052168
## Run 100 stress 0.1086814
## Run 101 stress 0.1052652
## Run 102 stress 0.1053421
## Run 103 stress 0.1055603
## Run 104 stress 0.1053419
## Run 105 stress 0.1049457
## Run 106 stress 0.1066368
## Run 107 stress 0.1064549
## Run 108 stress 0.1037063
## ... Procrustes: rmse 0.001999 max resid 0.02846922
## Run 109 stress 0.1051921
## Run 110 stress 0.1051284
## Run 111 stress 0.103972
## ... Procrustes: rmse 0.005570839 max resid 0.04323436
## Run 112 stress 0.1054927
## Run 113 stress 0.1054434
## Run 114 stress 0.1056562
## Run 115 stress 0.1066347
## Run 116 stress 0.1057859
## Run 117 stress 0.1051682
## Run 118 stress 0.1066037
## Run 119 stress 0.1053352
## Run 120 stress 0.1054651
## Run 121 stress 0.1055024
## Run 122 stress 0.1040743
## ... Procrustes: rmse 0.007981527 max resid 0.109647
## Run 123 stress 0.1053739
## Run 124 stress 0.1053455
## Run 125 stress 0.106802
## Run 126 stress 0.1039773
## ... Procrustes: rmse 0.005982936 max resid 0.0510808
## Run 127 stress 0.104097
## ... Procrustes: rmse 0.008873997 max resid 0.1236601
## Run 128 stress 0.1062863
## Run 129 stress 0.104212
## Run 130 stress 0.1053634
## Run 131 stress 0.1071667
## Run 132 stress 0.1037409
## ... Procrustes: rmse 0.002918317 max resid 0.02096524
## Run 133 stress 0.1051314
## Run 134 stress 0.1068572
## Run 135 stress 0.1051954
## Run 136 stress 0.1053888
## Run 137 stress 0.1069575
## Run 138 stress 0.1054316
## Run 139 stress 0.1082515
## Run 140 stress 0.1053717
## Run 141 stress 0.1065899
## Run 142 stress 0.1039
## ... Procrustes: rmse 0.005047809 max resid 0.06671632
## Run 143 stress 0.1037197
## ... Procrustes: rmse 0.002342273 max resid 0.02938035
## Run 144 stress 0.103986
## ... Procrustes: rmse 0.007920068 max resid 0.1183227
## Run 145 stress 0.1051985
## Run 146 stress 0.1039962
## ... Procrustes: rmse 0.008331081 max resid 0.124856
## Run 147 stress 0.1069593
## Run 148 stress 0.1039855
## ... Procrustes: rmse 0.008161781 max resid 0.1210844
## Run 149 stress 0.1055484
## Run 150 stress 0.1055342
## Run 151 stress 0.103879
## ... Procrustes: rmse 0.004546789 max resid 0.0431816
## Run 152 stress 0.1087563
## Run 153 stress 0.1072754
## Run 154 stress 0.105523
## Run 155 stress 0.1065537
## Run 156 stress 0.1068461
## Run 157 stress 0.1050111
## Run 158 stress 0.1050574
## Run 159 stress 0.1038874
## ... Procrustes: rmse 0.004471144 max resid 0.04615632
## Run 160 stress 0.1055009
## Run 161 stress 0.1054307
## Run 162 stress 0.1065876
## Run 163 stress 0.1040893
## ... Procrustes: rmse 0.008818907 max resid 0.1233156
## Run 164 stress 0.1055367
## Run 165 stress 0.1063708
## Run 166 stress 0.1038614
## ... Procrustes: rmse 0.005409715 max resid 0.07903686
## Run 167 stress 0.104871
## Run 168 stress 0.1055752
## Run 169 stress 0.1038039
## ... Procrustes: rmse 0.003589851 max resid 0.02210741
## Run 170 stress 0.1039836
## ... Procrustes: rmse 0.007939553 max resid 0.1172103
## Run 171 stress 0.1037232
## ... Procrustes: rmse 0.001996824 max resid 0.02116194
## Run 172 stress 0.1053043
## Run 173 stress 0.1051937
## Run 174 stress 0.1049951
## Run 175 stress 0.104861
## Run 176 stress 0.1036858
## ... Procrustes: rmse 0.0008393298 max resid 0.01193916
## Run 177 stress 0.1064913
## Run 178 stress 0.1038185
## ... Procrustes: rmse 0.003540973 max resid 0.0421968
## Run 179 stress 0.1044807
## Run 180 stress 0.103804
## ... Procrustes: rmse 0.003860474 max resid 0.0430969
## Run 181 stress 0.1037412
## ... Procrustes: rmse 0.003008108 max resid 0.02136756
## Run 182 stress 0.1066865
## Run 183 stress 0.1038053
## ... Procrustes: rmse 0.003558565 max resid 0.04302693
## Run 184 stress 0.1066373
## Run 185 stress 0.1048932
## Run 186 stress 0.1054485
## Run 187 stress 0.1037931
## ... Procrustes: rmse 0.003638536 max resid 0.04235056
## Run 188 stress 0.105176
## Run 189 stress 0.1049652
## Run 190 stress 0.103821
## ... Procrustes: rmse 0.004168878 max resid 0.04326228
## Run 191 stress 0.1039914
## ... Procrustes: rmse 0.006031939 max resid 0.07062848
## Run 192 stress 0.1051735
## Run 193 stress 0.1052762
## Run 194 stress 0.1052536
## Run 195 stress 0.1040649
## ... Procrustes: rmse 0.008499819 max resid 0.1191021
## Run 196 stress 0.1087645
## Run 197 stress 0.106407
## Run 198 stress 0.105535
## Run 199 stress 0.1086678
## Run 200 stress 0.1051626
## Run 201 stress 0.1052206
## Run 202 stress 0.1039534
## ... Procrustes: rmse 0.006070404 max resid 0.07777913
## Run 203 stress 0.1051645
## Run 204 stress 0.1057991
## Run 205 stress 0.1094136
## Run 206 stress 0.1054887
## Run 207 stress 0.1051517
## Run 208 stress 0.1051306
## Run 209 stress 0.1065712
## Run 210 stress 0.1040688
## ... Procrustes: rmse 0.008652522 max resid 0.1218283
## Run 211 stress 0.1078543
## Run 212 stress 0.1037077
## ... Procrustes: rmse 0.002066305 max resid 0.02862835
## Run 213 stress 0.1038549
## ... Procrustes: rmse 0.004129181 max resid 0.04233274
## Run 214 stress 0.1081731
## Run 215 stress 0.1054627
## Run 216 stress 0.1039943
## ... Procrustes: rmse 0.008283253 max resid 0.1240307
## Run 217 stress 0.1072255
## Run 218 stress 0.1053806
## Run 219 stress 0.1052059
## Run 220 stress 0.1051759
## Run 221 stress 0.1052599
## Run 222 stress 0.1102512
## Run 223 stress 0.1053471
## Run 224 stress 0.1078586
## Run 225 stress 0.105264
## Run 226 stress 0.1051973
## Run 227 stress 0.1042102
## Run 228 stress 0.1063842
## Run 229 stress 0.1040935
## ... Procrustes: rmse 0.007944401 max resid 0.1085664
## Run 230 stress 0.111074
## Run 231 stress 0.1069745
## Run 232 stress 0.1073528
## Run 233 stress 0.103821
## ... Procrustes: rmse 0.003117366 max resid 0.03421321
## Run 234 stress 0.1055506
## Run 235 stress 0.105613
## Run 236 stress 0.1039835
## ... Procrustes: rmse 0.007891517 max resid 0.1166033
## Run 237 stress 0.1049836
## Run 238 stress 0.1038503
## ... Procrustes: rmse 0.003974036 max resid 0.04246804
## Run 239 stress 0.1039039
## ... Procrustes: rmse 0.005193599 max resid 0.04323721
## Run 240 stress 0.1052643
## Run 241 stress 0.1056454
## Run 242 stress 0.1052456
## Run 243 stress 0.1051273
## Run 244 stress 0.103819
## ... Procrustes: rmse 0.004125337 max resid 0.04318105
## Run 245 stress 0.1039552
## ... Procrustes: rmse 0.005854683 max resid 0.06957821
## Run 246 stress 0.1054437
## Run 247 stress 0.1058083
## Run 248 stress 0.1053302
## Run 249 stress 0.1052671
## Run 250 stress 0.1070042
## Run 251 stress 0.1039844
## ... Procrustes: rmse 0.008057608 max resid 0.1196643
## Run 252 stress 0.1037111
## ... Procrustes: rmse 0.002050762 max resid 0.02656364
## Run 253 stress 0.1055361
## Run 254 stress 0.1040701
## ... Procrustes: rmse 0.008383744 max resid 0.116903
## Run 255 stress 0.105128
## Run 256 stress 0.1051727
## Run 257 stress 0.1058831
## Run 258 stress 0.105211
## Run 259 stress 0.1052793
## Run 260 stress 0.1038193
## ... Procrustes: rmse 0.004146237 max resid 0.0431907
## Run 261 stress 0.105438
## Run 262 stress 0.1050296
## Run 263 stress 0.1042594
## Run 264 stress 0.104895
## Run 265 stress 0.1101512
## Run 266 stress 0.1039903
## ... Procrustes: rmse 0.007480987 max resid 0.1095657
## Run 267 stress 0.103961
## ... Procrustes: rmse 0.005597371 max resid 0.04327574
## Run 268 stress 0.1049214
## Run 269 stress 0.1063384
## Run 270 stress 0.1039939
## ... Procrustes: rmse 0.00801165 max resid 0.1181229
## Run 271 stress 0.1060832
## Run 272 stress 0.1067616
## Run 273 stress 0.1051971
## Run 274 stress 0.1065053
## Run 275 stress 0.105126
## Run 276 stress 0.1083203
## Run 277 stress 0.105196
## Run 278 stress 0.1040693
## ... Procrustes: rmse 0.008306505 max resid 0.1179624
## Run 279 stress 0.1065383
## Run 280 stress 0.1055357
## Run 281 stress 0.1082019
## Run 282 stress 0.1055367
## Run 283 stress 0.1062985
## Run 284 stress 0.1038073
## ... Procrustes: rmse 0.003760729 max resid 0.05310769
## Run 285 stress 0.1037982
## ... Procrustes: rmse 0.003700509 max resid 0.04244458
## Run 286 stress 0.1040814
## ... Procrustes: rmse 0.008594948 max resid 0.1197492
## Run 287 stress 0.1052303
## Run 288 stress 0.103711
## ... Procrustes: rmse 0.00228947 max resid 0.03068486
## Run 289 stress 0.1065655
## Run 290 stress 0.1069532
## Run 291 stress 0.1056077
## Run 292 stress 0.1040067
## ... Procrustes: rmse 0.006302214 max resid 0.04339678
## Run 293 stress 0.1038044
## ... Procrustes: rmse 0.003888339 max resid 0.04310388
## Run 294 stress 0.106776
## Run 295 stress 0.1092293
## Run 296 stress 0.1039849
## ... Procrustes: rmse 0.008079487 max resid 0.1200314
## Run 297 stress 0.1053559
## Run 298 stress 0.1070089
## Run 299 stress 0.1041694
## ... Procrustes: rmse 0.009284797 max resid 0.1249693
## Run 300 stress 0.1054362
## *** Solution reached## 0.1
mds_whole_res_subset <- ps_nmds@sam_data %>%
as.tibble() %>%
select(Treatment,Temperature, Temperature_factor, Stressor,Stressor_numeric, Combinations, Temperature_factor, Lable, Name,F,G,M,A) %>%
bind_cols(as.tibble(scores(mds_whole_subset, display = "sites"))) %>%
mutate(T2=Temperature^2,T=Temperature_factor) %>%
select(-Temperature)## Warning in class(x) <- c(setdiff(subclass, tibble_class), tibble_class): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer be an
## S4 objectmds_whole_res_subset <- mds_whole_res_subset %>%
reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))
# NMDS subset anova combination
#NMDS1 anova combination
anova_NMDS1_mixed <- lm(NMDS1 ~F*G*M*A*T, data = mds_whole_res_subset)
autoplot(anova_NMDS1)
hist(residuals(anova_NMDS1 ))
anova(anova_NMDS1 )## Analysis of Variance Table
##
## Response: NMDS1
## Df Sum Sq Mean Sq F value Pr(>F)
## F 1 61.036 61.036 1028.3184 < 2.2e-16 ***
## G 1 0.002 0.002 0.0393 0.843076
## M 1 0.029 0.029 0.4889 0.485297
## A 1 2.625 2.625 44.2292 3.266e-10 ***
## T 2 0.652 0.326 5.4883 0.004844 **
## F:G 1 0.220 0.220 3.7031 0.055862 .
## F:M 1 0.000 0.000 0.0009 0.976463
## G:M 1 0.142 0.142 2.4002 0.123044
## F:A 1 3.009 3.009 50.7008 2.375e-11 ***
## G:A 1 0.014 0.014 0.2341 0.629091
## M:A 1 0.023 0.023 0.3956 0.530179
## F:T 2 2.659 1.329 22.3970 1.991e-09 ***
## G:T 2 0.068 0.034 0.5768 0.562723
## M:T 2 0.063 0.031 0.5270 0.591274
## A:T 2 0.211 0.105 1.7773 0.172005
## F:G:M 1 0.001 0.001 0.0203 0.886743
## F:G:A 1 0.128 0.128 2.1539 0.143927
## F:M:A 1 0.011 0.011 0.1884 0.664797
## G:M:A 1 0.130 0.130 2.1874 0.140864
## F:G:T 2 0.100 0.050 0.8402 0.433284
## F:M:T 2 0.033 0.017 0.2821 0.754527
## G:M:T 2 0.423 0.212 3.5659 0.030251 *
## F:A:T 2 0.452 0.226 3.8116 0.023888 *
## G:A:T 2 0.125 0.063 1.0551 0.350268
## M:A:T 2 0.037 0.019 0.3119 0.732470
## F:G:M:A 1 0.001 0.001 0.0243 0.876249
## F:G:M:T 2 0.046 0.023 0.3859 0.680386
## F:G:A:T 2 0.030 0.015 0.2526 0.777054
## F:M:A:T 2 0.016 0.008 0.1310 0.877296
## G:M:A:T 2 0.181 0.090 1.5243 0.220517
## F:G:M:A:T 2 0.021 0.011 0.1779 0.837146
## Residuals 183 10.862 0.059
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(anova_NMDS1)##
## Call:
## lm(formula = NMDS1 ~ F * G * M * A * T, data = mds_whole_res)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.90640 -0.12562 -0.00319 0.08877 0.98663
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.62357 0.10895 -5.723 4.21e-08 ***
## F 1.01436 0.15408 6.583 4.72e-10 ***
## G 0.08727 0.16343 0.534 0.5940
## M -0.10535 0.15408 -0.684 0.4950
## A 0.25032 0.15408 1.625 0.1060
## T24 0.04710 0.15408 0.306 0.7602
## T28 -0.15322 0.15408 -0.994 0.3213
## F:G -0.03967 0.22462 -0.177 0.8600
## F:M 0.15503 0.21791 0.711 0.4777
## G:M -0.07758 0.22462 -0.345 0.7302
## F:A -0.34377 0.22462 -1.531 0.1276
## G:A -0.08645 0.22462 -0.385 0.7008
## M:A 0.07985 0.21791 0.366 0.7145
## F:T24 0.11684 0.22462 0.520 0.6036
## F:T28 0.51957 0.21791 2.384 0.0181 *
## G:T24 -0.40315 0.22462 -1.795 0.0743 .
## G:T28 -0.08559 0.22462 -0.381 0.7036
## M:T24 -0.11534 0.21791 -0.529 0.5972
## M:T28 -0.02001 0.21791 -0.092 0.9269
## A:T24 0.30010 0.21791 1.377 0.1701
## A:T28 0.20718 0.21791 0.951 0.3430
## F:G:M -0.03538 0.31295 -0.113 0.9101
## F:G:A 0.32514 0.31765 1.024 0.3074
## F:M:A -0.12032 0.31295 -0.384 0.7011
## G:M:A 0.01457 0.31295 0.047 0.9629
## F:G:T24 0.19467 0.31765 0.613 0.5407
## F:G:T28 -0.01761 0.31295 -0.056 0.9552
## F:M:T24 -0.26511 0.31295 -0.847 0.3980
## F:M:T28 -0.11238 0.30817 -0.365 0.7158
## G:M:T24 0.62782 0.32535 1.930 0.0552 .
## G:M:T28 0.16828 0.31295 0.538 0.5914
## F:A:T24 -0.36116 0.32229 -1.121 0.2639
## F:A:T28 -0.20663 0.31295 -0.660 0.5099
## G:A:T24 0.18880 0.31765 0.594 0.5530
## G:A:T28 -0.05402 0.31295 -0.173 0.8631
## M:A:T24 0.20496 0.30817 0.665 0.5068
## M:A:T28 0.02853 0.30817 0.093 0.9263
## F:G:M:A -0.07491 0.44257 -0.169 0.8658
## F:G:M:T24 0.09541 0.45470 0.210 0.8340
## F:G:M:T28 0.03307 0.43921 0.075 0.9401
## F:G:A:T24 -0.33135 0.45252 -0.732 0.4650
## F:G:A:T28 -0.01474 0.44257 -0.033 0.9735
## F:M:A:T24 0.06044 0.44592 0.136 0.8923
## F:M:A:T28 0.19036 0.43921 0.433 0.6652
## G:M:A:T24 -0.62681 0.45143 -1.388 0.1667
## G:M:A:T28 0.02288 0.43921 0.052 0.9585
## F:G:M:A:T24 0.24337 0.63919 0.381 0.7038
## F:G:M:A:T28 -0.13143 0.62114 -0.212 0.8327
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2436 on 183 degrees of freedom
## Multiple R-squared: 0.8697, Adjusted R-squared: 0.8362
## F-statistic: 25.98 on 47 and 183 DF, p-value: < 2.2e-16#NMDS2 anova combination
anova_NMDS2_mixed <- lm(NMDS2 ~F*G*M*A*T, data = mds_whole_res_subset)
autoplot(anova_NMDS2)
hist(residuals(anova_NMDS2 ))
anova(anova_NMDS2 )## Analysis of Variance Table
##
## Response: NMDS2
## Df Sum Sq Mean Sq F value Pr(>F)
## F 1 0.0493 0.0493 0.3960 0.529963
## G 1 0.0004 0.0004 0.0034 0.953514
## M 1 0.7552 0.7552 6.0673 0.014695 *
## A 1 4.8989 4.8989 39.3559 2.484e-09 ***
## T 2 7.7550 3.8775 31.1500 2.275e-12 ***
## F:G 1 0.1971 0.1971 1.5831 0.209918
## F:M 1 0.0523 0.0523 0.4200 0.517749
## G:M 1 0.0059 0.0059 0.0471 0.828501
## F:A 1 0.9378 0.9378 7.5338 0.006659 **
## G:A 1 0.4115 0.4115 3.3061 0.070658 .
## M:A 1 0.1714 0.1714 1.3766 0.242201
## F:T 2 0.0725 0.0362 0.2911 0.747764
## G:T 2 0.2108 0.1054 0.8467 0.430497
## M:T 2 0.2787 0.1394 1.1196 0.328627
## A:T 2 0.2865 0.1433 1.1508 0.318652
## F:G:M 1 0.0888 0.0888 0.7132 0.399501
## F:G:A 1 0.0137 0.0137 0.1101 0.740426
## F:M:A 1 0.3646 0.3646 2.9292 0.088686 .
## G:M:A 1 0.4031 0.4031 3.2382 0.073587 .
## F:G:T 2 0.0916 0.0458 0.3678 0.692757
## F:M:T 2 0.1976 0.0988 0.7937 0.453731
## G:M:T 2 0.3708 0.1854 1.4896 0.228179
## F:A:T 2 1.1742 0.5871 4.7163 0.010064 *
## G:A:T 2 0.4435 0.2217 1.7813 0.171333
## M:A:T 2 0.1121 0.0561 0.4503 0.638123
## F:G:M:A 1 0.0078 0.0078 0.0624 0.803066
## F:G:M:T 2 0.0597 0.0299 0.2399 0.786919
## F:G:A:T 2 0.2561 0.1281 1.0287 0.359528
## F:M:A:T 2 0.0270 0.0135 0.1084 0.897316
## G:M:A:T 2 0.1251 0.0625 0.5025 0.605855
## F:G:M:A:T 2 0.7768 0.3884 3.1202 0.046506 *
## Residuals 183 22.7795 0.1245
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(anova_NMDS2)##
## Call:
## lm(formula = NMDS2 ~ F * G * M * A * T, data = mds_whole_res)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.67576 -0.10981 -0.00263 0.11760 1.33865
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.159311 0.157783 -1.010 0.3140
## F 0.085881 0.223139 0.385 0.7008
## G 0.067201 0.236675 0.284 0.7768
## M 0.126427 0.223139 0.567 0.5717
## A -0.231951 0.223139 -1.039 0.2999
## T24 0.434002 0.223139 1.945 0.0533 .
## T28 0.554193 0.223139 2.484 0.0139 *
## F:G 0.393576 0.325279 1.210 0.2279
## F:M -0.274407 0.315567 -0.870 0.3857
## G:M -0.153747 0.325279 -0.473 0.6370
## F:A -0.038816 0.325279 -0.119 0.9051
## G:A 0.076456 0.325279 0.235 0.8144
## M:A -0.027077 0.315567 -0.086 0.9317
## F:T24 -0.103381 0.325279 -0.318 0.7510
## F:T28 -0.230081 0.315567 -0.729 0.4669
## G:T24 0.144518 0.325279 0.444 0.6574
## G:T28 -0.032296 0.325279 -0.099 0.9210
## M:T24 -0.144553 0.315567 -0.458 0.6474
## M:T28 -0.088370 0.315567 -0.280 0.7798
## A:T24 0.094618 0.315567 0.300 0.7646
## A:T28 -0.167913 0.315567 -0.532 0.5953
## F:G:M -0.335092 0.453198 -0.739 0.4606
## F:G:A -0.647511 0.460014 -1.408 0.1609
## F:M:A 0.149664 0.453198 0.330 0.7416
## G:M:A -0.017552 0.453198 -0.039 0.9691
## F:G:T24 -0.742368 0.460014 -1.614 0.1083
## F:G:T28 -0.181871 0.453198 -0.401 0.6887
## F:M:T24 -0.227877 0.453198 -0.503 0.6157
## F:M:T28 0.192243 0.446279 0.431 0.6671
## G:M:T24 -0.097784 0.471154 -0.208 0.8358
## G:M:T28 0.031290 0.453198 0.069 0.9450
## F:A:T24 -0.009746 0.466730 -0.021 0.9834
## F:A:T28 0.541322 0.453198 1.194 0.2338
## G:A:T24 -1.011734 0.460014 -2.199 0.0291 *
## G:A:T28 -0.090973 0.453198 -0.201 0.8411
## M:A:T24 -0.449489 0.446279 -1.007 0.3152
## M:A:T28 -0.044948 0.446279 -0.101 0.9199
## F:G:M:A 0.580079 0.640919 0.905 0.3666
## F:G:M:T24 1.183700 0.658482 1.798 0.0739 .
## F:G:M:T28 0.151148 0.636045 0.238 0.8124
## F:G:A:T24 1.576498 0.655324 2.406 0.0171 *
## F:G:A:T28 0.138070 0.640919 0.215 0.8297
## F:M:A:T24 0.691480 0.645756 1.071 0.2857
## F:M:A:T28 -0.283331 0.636045 -0.445 0.6565
## G:M:A:T24 0.766263 0.653739 1.172 0.2427
## G:M:A:T28 0.153726 0.636045 0.242 0.8093
## F:G:M:A:T24 -1.853952 0.925648 -2.003 0.0467 *
## F:G:M:A:T28 0.279283 0.899504 0.310 0.7565
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3528 on 183 degrees of freedom
## Multiple R-squared: 0.4748, Adjusted R-squared: 0.3399
## F-statistic: 3.52 on 47 and 183 DF, p-value: 6.71e-10#NMDS3 anova combination
anova_NMDS3 <- lm(NMDS3 ~F*G*M*A*T, data = mds_whole_res_subset)
autoplot(anova_NMDS3)
hist(residuals(anova_NMDS3 ))
anova(anova_NMDS3 )## Analysis of Variance Table
##
## Response: NMDS3
## Df Sum Sq Mean Sq F value Pr(>F)
## F 1 0.1660 0.16603 1.1221 0.29087
## G 1 0.0096 0.00956 0.0646 0.79969
## M 1 0.0000 0.00005 0.0003 0.98601
## A 1 0.2031 0.20307 1.3724 0.24292
## T 2 3.9849 1.99245 13.4656 3.501e-06 ***
## F:G 1 0.0296 0.02957 0.1998 0.65537
## F:M 1 0.4154 0.41544 2.8076 0.09552 .
## G:M 1 0.4513 0.45134 3.0503 0.08240 .
## F:A 1 0.2174 0.21738 1.4691 0.22705
## G:A 1 0.9208 0.92075 6.2227 0.01350 *
## M:A 1 0.3730 0.37299 2.5208 0.11408
## F:T 2 5.6642 2.83212 19.1404 2.831e-08 ***
## G:T 2 0.0494 0.02469 0.1669 0.84645
## M:T 2 0.4793 0.23963 1.6195 0.20082
## A:T 2 0.3679 0.18394 1.2432 0.29090
## F:G:M 1 0.4262 0.42623 2.8806 0.09135 .
## F:G:A 1 0.0048 0.00481 0.0325 0.85716
## F:M:A 1 0.0165 0.01645 0.1112 0.73917
## G:M:A 1 0.2820 0.28196 1.9056 0.16914
## F:G:T 2 0.3142 0.15709 1.0617 0.34799
## F:M:T 2 0.0512 0.02560 0.1730 0.84125
## G:M:T 2 0.0206 0.01028 0.0695 0.93293
## F:A:T 2 0.0221 0.01107 0.0748 0.92793
## G:A:T 2 0.1762 0.08810 0.5954 0.55242
## M:A:T 2 0.4557 0.22784 1.5398 0.21719
## F:G:M:A 1 0.1968 0.19683 1.3302 0.25027
## F:G:M:T 2 0.2757 0.13783 0.9315 0.39583
## F:G:A:T 2 0.0409 0.02044 0.1382 0.87105
## F:M:A:T 2 0.0524 0.02619 0.1770 0.83792
## G:M:A:T 2 1.2895 0.64477 4.3575 0.01417 *
## F:G:M:A:T 2 0.4549 0.22747 1.5373 0.21772
## Residuals 183 27.0778 0.14797
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(anova_NMDS3)##
## Call:
## lm(formula = NMDS3 ~ F * G * M * A * T, data = mds_whole_res_subset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.08527 -0.18913 -0.00264 0.14258 1.25117
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.42168 0.17203 2.451 0.01517 *
## F -0.43074 0.24328 -1.771 0.07830 .
## G -0.05715 0.25804 -0.221 0.82498
## M -0.35646 0.24328 -1.465 0.14458
## A -0.07423 0.24328 -0.305 0.76063
## T24 -0.50528 0.24328 -2.077 0.03921 *
## T28 -0.77280 0.24328 -3.177 0.00175 **
## F:G 0.32615 0.35464 0.920 0.35896
## F:M 0.15911 0.34405 0.462 0.64429
## G:M 0.75950 0.35464 2.142 0.03355 *
## F:A -0.02962 0.35464 -0.084 0.93354
## G:A -0.06466 0.35464 -0.182 0.85553
## M:A 0.61925 0.34405 1.800 0.07353 .
## F:T24 0.61182 0.35464 1.725 0.08618 .
## F:T28 0.48326 0.34405 1.405 0.16183
## G:T24 -0.10935 0.35464 -0.308 0.75817
## G:T28 0.09818 0.35464 0.277 0.78222
## M:T24 0.21892 0.34405 0.636 0.52538
## M:T28 0.38624 0.34405 1.123 0.26307
## A:T24 0.40884 0.34405 1.188 0.23625
## A:T28 0.31402 0.34405 0.913 0.36260
## F:G:M -0.99117 0.49411 -2.006 0.04633 *
## F:G:A -0.36071 0.50154 -0.719 0.47293
## F:M:A -0.21900 0.49411 -0.443 0.65813
## G:M:A -0.69423 0.49411 -1.405 0.16171
## F:G:T24 -0.02092 0.50154 -0.042 0.96677
## F:G:T28 0.09487 0.49411 0.192 0.84796
## F:M:T24 -0.24707 0.49411 -0.500 0.61765
## F:M:T28 0.03690 0.48657 0.076 0.93962
## G:M:T24 -0.47948 0.51369 -0.933 0.35183
## G:M:T28 -0.84546 0.49411 -1.711 0.08876 .
## F:A:T24 -0.34859 0.50886 -0.685 0.49419
## F:A:T28 0.34959 0.49411 0.708 0.48015
## G:A:T24 -0.27441 0.50154 -0.547 0.58496
## G:A:T28 -0.27332 0.49411 -0.553 0.58084
## M:A:T24 -0.86023 0.48657 -1.768 0.07873 .
## M:A:T28 -0.63003 0.48657 -1.295 0.19700
## F:G:M:A 0.83301 0.69878 1.192 0.23476
## F:G:M:T24 0.76706 0.71792 1.068 0.28673
## F:G:M:T28 0.50833 0.69346 0.733 0.46448
## F:G:A:T24 0.65627 0.71448 0.919 0.35955
## F:G:A:T28 -0.35146 0.69878 -0.503 0.61559
## F:M:A:T24 0.59664 0.70405 0.847 0.39785
## F:M:A:T28 -0.41183 0.69346 -0.594 0.55333
## G:M:A:T24 0.96117 0.71275 1.349 0.17916
## G:M:A:T28 1.25184 0.69346 1.805 0.07269 .
## F:G:M:A:T24 -1.41279 1.00921 -1.400 0.16323
## F:G:M:A:T28 0.22341 0.98070 0.228 0.82005
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3847 on 183 degrees of freedom
## Multiple R-squared: 0.3914, Adjusted R-squared: 0.235
## F-statistic: 2.504 on 47 and 183 DF, p-value: 7.316e-06#NMDS subset anova numeric
#NMDS3
anova_NMDS3_numeric <- lm(NMDS3 ~ Stressor_numeric*T, data = mds_whole_res_subset)
autoplot(anova_NMDS3_numeric)
hist(residuals(anova_NMDS3_numeric ))
anova(anova_NMDS3_numeric )## Analysis of Variance Table
##
## Response: NMDS3
## Df Sum Sq Mean Sq F value Pr(>F)
## Stressor_numeric 1 0.004 0.00405 0.0236 0.878077
## T 2 3.979 1.98941 11.5879 1.622e-05 ***
## Stressor_numeric:T 2 1.878 0.93918 5.4705 0.004788 **
## Residuals 225 38.628 0.17168
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(anova_NMDS3_numeric)##
## Call:
## lm(formula = NMDS3 ~ Stressor_numeric * T, data = mds_whole_res_subset)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17425 -0.22735 -0.04243 0.17112 1.46893
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.38990 0.10491 3.716 0.000255 ***
## Stressor_numeric -0.10677 0.04662 -2.290 0.022936 *
## T24 -0.43546 0.14964 -2.910 0.003977 **
## T28 -0.74761 0.14743 -5.071 8.27e-07 ***
## Stressor_numeric:T24 0.11091 0.06705 1.654 0.099482 .
## Stressor_numeric:T28 0.21739 0.06572 3.308 0.001095 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4143 on 225 degrees of freedom
## Multiple R-squared: 0.1317, Adjusted R-squared: 0.1125
## F-statistic: 6.828 on 5 and 225 DF, p-value: 5.956e-06#means of NMDS
means_NMDS_subset <- mds_whole_res_subset %>%
filter(Temperature_factor %in% c("20", "24", "28")) %>%
dplyr::group_by(Treatment,Combinations, Stressor, Stressor_numeric, Temperature_factor,F,G,M,A) %>%
dplyr::summarize(NMDS1_mean = mean(NMDS1),
NMDS2_mean = mean(NMDS2),
NMDS3_mean = mean(NMDS3),
NMDS1_SE = stderr(NMDS1),
NMDS2_SE = stderr(NMDS2),
NMDS3_SE= stderr(NMDS3))## `summarise()` has grouped output by 'Treatment', 'Combinations', 'Stressor', 'Stressor_numeric', 'Temperature_factor', 'F', 'G', 'M'. You can override using the `.groups` argument.#NMDS1 subpopulation
nmds1_plot_temperature_subset <- means_NMDS_subset %>%
ggplot() +
geom_point(aes(x=Temperature_factor, y=NMDS1_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Temperature_factor,y=NMDS1_mean,col=Combinations,ymin=NMDS1_mean-NMDS1_SE, ymax=NMDS1_mean+NMDS1_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
geom_line(aes(group=Combinations, x=Temperature_factor, y=NMDS1_mean, col=Combinations),position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
scale_colour_viridis_d(direction = -1)+
xlab("Temperature") +
ylab("NMDS1_mean subpopulation")+
theme(legend.position = "none") +
ylab("NMDS1_subpopulation (mean)")
plot_NMDS1_means_numeric_subset<- means_NMDS_subset %>%
ggplot() +
geom_point(aes(x=Stressor_numeric, y=NMDS1_mean, col=Temperature_factor, shape= Temperature_factor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=NMDS1_mean,col=Temperature_factor,ymin=NMDS1_mean-NMDS1_SE, ymax=NMDS1_mean+NMDS1_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
ylab("NMDS1_subpopulation (mean)") +
xlab( "Number of stressors")+
theme(legend.position = "none")
# plot
nmds1_plot_temperature_subset + plot_NMDS1_means_numeric_subset 
#NMDS2 subpopulation
NMDS2_plot_temperature_subset <- means_NMDS_subset %>%
ggplot() +
geom_point(aes(x=Temperature_factor, y=NMDS2_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Temperature_factor,y=NMDS2_mean,col=Combinations,ymin=NMDS2_mean-NMDS2_SE, ymax=NMDS2_mean+NMDS2_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
geom_line(aes(group=Combinations, x=Temperature_factor, y=NMDS2_mean, col=Combinations),position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
scale_colour_viridis_d(direction = -1)+
xlab("Temperature") +
theme(legend.position = "none") +
ylab("NMDS2_subpopulation (mean)")
plot_NMDS2_means_numeric_subset <- means_NMDS_subset %>%
ggplot() +
geom_point(aes(x=Stressor_numeric, y=NMDS2_mean, col=Temperature_factor, shape= Temperature_factor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=NMDS2_mean,col=Temperature_factor,ymin=NMDS2_mean-NMDS2_SE, ymax=NMDS2_mean+NMDS2_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
ylab("NMDS2_subpopulation (mean)") +
xlab( "Number of stressors")+
theme(legend.position = "none")
#plot
NMDS2_plot_temperature_subset + plot_NMDS2_means_numeric_subset 
Dataset of single genera, namely Tychonema_CCAP_1459-11B, Sulfuricurvum and Chromatium, calculating NMDS for all, respectively
ps_relative_cyano <- subset_taxa(ps_relative, Genus ==
"Tychonema_CCAP_1459-11B")
df_cyano <- psmelt(ps_relative_cyano)
df_cyano <- df_cyano %>%
mutate(T2=Temperature^2,T=Temperature_factor) %>%
select(-Temperature)
df_cyano <- df_cyano %>%
filter(Treatment != "ori")
df_cyano <- df_cyano %>%
dplyr::group_by(Treatment,Stressor_numeric,Lable,Combinations, T,F,G,M,A) %>%
dplyr::summarize(Abundance = sum(Abundance)) %>%
reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))## `summarise()` has grouped output by 'Treatment', 'Stressor_numeric', 'Lable', 'Combinations', 'T', 'F', 'G', 'M'. You can override using the `.groups` argument.#anova on combination
anova_tychonea <- lm(logit(Abundance) ~F*G*M*A*T, data = df_cyano)
autoplot(anova_tychonea)
hist(residuals(anova_tychonea))
anova(anova_tychonea)## Analysis of Variance Table
##
## Response: logit(Abundance)
## Df Sum Sq Mean Sq F value Pr(>F)
## F 1 93.372 93.372 191.0631 < 2.2e-16 ***
## G 1 0.037 0.037 0.0750 0.784453
## M 1 0.304 0.304 0.6214 0.431543
## A 1 79.227 79.227 162.1187 < 2.2e-16 ***
## T 2 4.268 2.134 4.3663 0.014047 *
## F:G 1 0.062 0.062 0.1274 0.721557
## F:M 1 0.014 0.014 0.0277 0.868010
## G:M 1 2.988 2.988 6.1141 0.014323 *
## F:A 1 77.297 77.297 158.1699 < 2.2e-16 ***
## G:A 1 0.391 0.391 0.8007 0.372060
## M:A 1 2.395 2.395 4.9007 0.028081 *
## F:T 2 6.685 3.342 6.8392 0.001366 **
## G:T 2 2.332 1.166 2.3864 0.094817 .
## M:T 2 2.128 1.064 2.1775 0.116249
## A:T 2 1.412 0.706 1.4450 0.238412
## F:G:M 1 0.655 0.655 1.3403 0.248492
## F:G:A 1 0.042 0.042 0.0859 0.769814
## F:M:A 1 0.042 0.042 0.0851 0.770875
## G:M:A 1 0.037 0.037 0.0767 0.782162
## F:G:T 2 2.220 1.110 2.2709 0.106125
## F:M:T 2 0.728 0.364 0.7443 0.476484
## G:M:T 2 2.575 1.288 2.6346 0.074468 .
## F:A:T 2 0.104 0.052 0.1067 0.898876
## G:A:T 2 0.305 0.152 0.3120 0.732335
## M:A:T 2 0.976 0.488 0.9985 0.370426
## F:G:M:A 1 0.474 0.474 0.9696 0.326080
## F:G:M:T 2 0.139 0.070 0.1423 0.867456
## F:G:A:T 2 0.202 0.101 0.2068 0.813350
## F:M:A:T 2 0.338 0.169 0.3454 0.708425
## G:M:A:T 2 1.333 0.667 1.3638 0.258261
## F:G:M:A:T 2 0.842 0.421 0.8612 0.424347
## Residuals 183 89.432 0.489
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(anova_tychonea)##
## Call:
## lm(formula = logit(Abundance) ~ F * G * M * A * T, data = df_cyano)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.48855 -0.31281 -0.00915 0.29578 2.41578
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.10233 0.31263 -6.725 2.18e-10 ***
## F -2.08698 0.44213 -4.720 4.67e-06 ***
## G -0.38890 0.46895 -0.829 0.408013
## M 0.59739 0.44213 1.351 0.178310
## A -2.19776 0.44213 -4.971 1.53e-06 ***
## T24 -0.57596 0.44213 -1.303 0.194319
## T28 0.13664 0.44213 0.309 0.757635
## F:G 0.04950 0.64451 0.077 0.938860
## F:M -0.52908 0.62527 -0.846 0.398566
## G:M -0.76128 0.64451 -1.181 0.239068
## F:A 2.16353 0.64451 3.357 0.000959 ***
## G:A -0.12884 0.64451 -0.200 0.841776
## M:A -1.05617 0.62527 -1.689 0.092893 .
## F:T24 0.11824 0.64451 0.183 0.854638
## F:T28 -0.41784 0.62527 -0.668 0.504815
## G:T24 1.65914 0.64451 2.574 0.010836 *
## G:T28 0.60803 0.64451 0.943 0.346723
## M:T24 0.84276 0.62527 1.348 0.179378
## M:T28 -0.26699 0.62527 -0.427 0.669885
## A:T24 0.36839 0.62527 0.589 0.556470
## A:T28 0.01890 0.62527 0.030 0.975913
## F:G:M 1.19959 0.89797 1.336 0.183243
## F:G:A 0.57908 0.91147 0.635 0.526009
## F:M:A 0.78792 0.89797 0.877 0.381396
## G:M:A 1.35218 0.89797 1.506 0.133837
## F:G:T24 -0.91244 0.91147 -1.001 0.318119
## F:G:T28 -0.50304 0.89797 -0.560 0.576027
## F:M:T24 -0.28087 0.89797 -0.313 0.754805
## F:M:T28 0.38004 0.88426 0.430 0.667859
## G:M:T24 -1.09066 0.93355 -1.168 0.244208
## G:M:T28 0.75242 0.89797 0.838 0.403174
## F:A:T24 0.01639 0.92478 0.018 0.985883
## F:A:T28 -0.08825 0.89797 -0.098 0.921815
## G:A:T24 -0.40545 0.91147 -0.445 0.656968
## G:A:T28 0.72491 0.89797 0.807 0.420557
## M:A:T24 0.51831 0.88426 0.586 0.558497
## M:A:T28 0.53227 0.88426 0.602 0.547959
## F:G:M:A -1.94803 1.26992 -1.534 0.126762
## F:G:M:T24 -0.19651 1.30472 -0.151 0.880446
## F:G:M:T28 -1.05818 1.26026 -0.840 0.402200
## F:G:A:T24 -0.11091 1.29846 -0.085 0.932025
## F:G:A:T28 -0.85559 1.26992 -0.674 0.501332
## F:M:A:T24 -0.50726 1.27951 -0.396 0.692235
## F:M:A:T28 -0.45012 1.26026 -0.357 0.721381
## G:M:A:T24 -0.32512 1.29532 -0.251 0.802097
## G:M:A:T28 -2.25891 1.26026 -1.792 0.074719 .
## F:G:M:A:T24 1.31007 1.83408 0.714 0.475958
## F:G:M:A:T28 2.33497 1.78228 1.310 0.191804
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6991 on 183 degrees of freedom
## Multiple R-squared: 0.7605, Adjusted R-squared: 0.6989
## F-statistic: 12.36 on 47 and 183 DF, p-value: < 2.2e-16# Sulfuricurvum
ps_sulfuri <- subset_taxa(ps_relative, Genus == "Sulfuricurvum")
df_sulfuri <- psmelt(ps_sulfuri)
df_sulfuri <- df_sulfuri %>%
mutate(T2=Temperature^2,T=Temperature_factor) %>%
select(-Temperature)
df_sulfuri <- df_sulfuri %>%
filter(Treatment!="ori") %>%
dplyr::group_by(Treatment,Stressor_numeric,Combinations, Lable,T,F,G,M,A) %>%
dplyr::summarize(Abundance = sum(Abundance))## `summarise()` has grouped output by 'Treatment', 'Stressor_numeric', 'Combinations', 'Lable', 'T', 'F', 'G', 'M'. You can override using the `.groups` argument.#anova on combinations
anova_sulfuri <- lm(asin(sqrt(Abundance)) ~F*G*M*A*T, data = df_sulfuri)
autoplot(anova_sulfuri)
hist(residuals(anova_sulfuri))
anova(anova_sulfuri)## Analysis of Variance Table
##
## Response: asin(sqrt(Abundance))
## Df Sum Sq Mean Sq F value Pr(>F)
## F 1 0.03216 0.032159 11.0700 0.0010604 **
## G 1 0.04274 0.042741 14.7128 0.0001722 ***
## M 1 0.00026 0.000256 0.0880 0.7670137
## A 1 0.03900 0.038996 13.4236 0.0003253 ***
## T 2 0.03474 0.017369 5.9790 0.0030527 **
## F:G 1 0.00654 0.006539 2.2509 0.1352607
## F:M 1 0.00647 0.006475 2.2287 0.1371887
## G:M 1 0.00015 0.000153 0.0527 0.8187109
## F:A 1 0.06559 0.065589 22.5777 4.069e-06 ***
## G:A 1 0.00159 0.001586 0.5460 0.4608872
## M:A 1 0.00258 0.002580 0.8883 0.3471903
## F:T 2 0.02910 0.014550 5.0084 0.0076261 **
## G:T 2 0.00094 0.000471 0.1621 0.8504703
## M:T 2 0.01638 0.008189 2.8188 0.0622686 .
## A:T 2 0.03860 0.019301 6.6439 0.0016388 **
## F:G:M 1 0.00183 0.001828 0.6293 0.4286264
## F:G:A 1 0.00713 0.007128 2.4535 0.1189889
## F:M:A 1 0.00007 0.000071 0.0243 0.8763079
## G:M:A 1 0.00813 0.008128 2.7978 0.0961047 .
## F:G:T 2 0.00806 0.004032 1.3878 0.2522327
## F:M:T 2 0.02019 0.010097 3.4757 0.0329978 *
## G:M:T 2 0.00117 0.000583 0.2007 0.8183574
## F:A:T 2 0.00100 0.000499 0.1718 0.8423117
## G:A:T 2 0.00371 0.001857 0.6392 0.5289039
## M:A:T 2 0.01492 0.007458 2.5673 0.0795012 .
## F:G:M:A 1 0.00034 0.000339 0.1166 0.7331165
## F:G:M:T 2 0.01309 0.006547 2.2538 0.1079070
## F:G:A:T 2 0.00449 0.002247 0.7735 0.4629033
## F:M:A:T 2 0.00230 0.001151 0.3961 0.6735055
## G:M:A:T 2 0.00395 0.001973 0.6791 0.5083566
## F:G:M:A:T 2 0.00374 0.001871 0.6440 0.5263906
## Residuals 183 0.53162 0.002905
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(anova_sulfuri)##
## Call:
## lm(formula = asin(sqrt(Abundance)) ~ F * G * M * A * T, data = df_sulfuri)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.144494 -0.020506 -0.005465 0.015081 0.259963
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.0708255 0.0241041 2.938 0.00372 **
## F -0.0067231 0.0340884 -0.197 0.84387
## G 0.0074519 0.0361562 0.206 0.83694
## M -0.0294875 0.0340884 -0.865 0.38815
## A -0.0503197 0.0340884 -1.476 0.14162
## T24 -0.0543332 0.0340884 -1.594 0.11269
## T28 -0.0654043 0.0340884 -1.919 0.05658 .
## F:G 0.0654081 0.0496919 1.316 0.18973
## F:M 0.0608374 0.0482083 1.262 0.20857
## G:M 0.0000473 0.0496919 0.001 0.99924
## F:A -0.0055650 0.0496919 -0.112 0.91095
## G:A 0.0255688 0.0496919 0.515 0.60749
## M:A 0.0384315 0.0482083 0.797 0.42637
## F:T24 0.1123383 0.0496919 2.261 0.02496 *
## F:T28 0.0187107 0.0482083 0.388 0.69838
## G:T24 0.0041848 0.0496919 0.084 0.93298
## G:T28 0.0224243 0.0496919 0.451 0.65233
## M:T24 0.0524766 0.0482083 1.089 0.27779
## M:T28 0.0495668 0.0482083 1.028 0.30522
## A:T24 0.0956222 0.0482083 1.984 0.04880 *
## A:T28 0.0620968 0.0482083 1.288 0.19934
## F:G:M -0.0321723 0.0692338 -0.465 0.64271
## F:G:A -0.1011814 0.0702750 -1.440 0.15163
## F:M:A -0.0779992 0.0692338 -1.127 0.26138
## G:M:A -0.0365612 0.0692338 -0.528 0.59808
## F:G:T24 -0.0322423 0.0702750 -0.459 0.64692
## F:G:T28 0.0118470 0.0692338 0.171 0.86432
## F:M:T24 -0.1272817 0.0692338 -1.838 0.06762 .
## F:M:T28 -0.0604214 0.0681767 -0.886 0.37665
## G:M:T24 -0.0178443 0.0719768 -0.248 0.80448
## G:M:T28 -0.0195382 0.0692338 -0.282 0.77810
## F:A:T24 -0.0916130 0.0713010 -1.285 0.20046
## F:A:T28 -0.0091811 0.0692338 -0.133 0.89465
## G:A:T24 -0.0301866 0.0702750 -0.430 0.66803
## G:A:T28 -0.0631096 0.0692338 -0.912 0.36321
## M:A:T24 -0.0900777 0.0681767 -1.321 0.18807
## M:A:T28 -0.0278244 0.0681767 -0.408 0.68366
## F:G:M:A 0.1082929 0.0979114 1.106 0.27017
## F:G:M:T24 0.0555469 0.1005944 0.552 0.58149
## F:G:M:T28 -0.0500372 0.0971668 -0.515 0.60720
## F:G:A:T24 0.1073489 0.1001119 1.072 0.28500
## F:G:A:T28 0.0362914 0.0979114 0.371 0.71132
## F:M:A:T24 0.1125290 0.0986504 1.141 0.25549
## F:M:A:T28 0.0792571 0.0971668 0.816 0.41574
## G:M:A:T24 0.0671390 0.0998698 0.672 0.50226
## G:M:A:T28 0.1496609 0.0971668 1.540 0.12523
## F:G:M:A:T24 -0.1123707 0.1414086 -0.795 0.42784
## F:G:M:A:T28 -0.1507164 0.1374146 -1.097 0.27417
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0539 on 183 degrees of freedom
## Multiple R-squared: 0.436, Adjusted R-squared: 0.2911
## F-statistic: 3.01 on 47 and 183 DF, p-value: 7.15e-08Abiotic factors and oxygen
abiotic<- ps_relative@sam_data %>%
as.tibble() %>%
select(Treatment,Temperature, Combinations, Temperature_factor, Lable, Name,F,G,M,A, Oxygen, TN,TC, pH,Stressor_numeric, Stressor) %>%
mutate(T2=Temperature^2,T=Temperature_factor) %>%
select(-Temperature)## Warning in class(x) <- c(setdiff(subclass, tibble_class), tibble_class): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer be an
## S4 objectpH_anova <- lm(pH ~F*G*M*A*T, data = abiotic)
autoplot(pH_anova)
hist(residuals(pH_anova))
anova(pH_anova)## Analysis of Variance Table
##
## Response: pH
## Df Sum Sq Mean Sq F value Pr(>F)
## F 1 357.64 357.64 3216.0251 < 2.2e-16 ***
## G 1 0.14 0.14 1.2477 0.2654646
## M 1 0.10 0.10 0.8570 0.3557852
## A 1 1.77 1.77 15.9178 9.559e-05 ***
## T 2 0.83 0.42 3.7372 0.0256572 *
## F:G 1 0.03 0.03 0.2885 0.5918137
## F:M 1 0.27 0.27 2.3868 0.1240949
## G:M 1 0.07 0.07 0.6723 0.4133055
## F:A 1 1.62 1.62 14.6075 0.0001813 ***
## G:A 1 0.15 0.15 1.3826 0.2411803
## M:A 1 0.43 0.43 3.9011 0.0497571 *
## F:T 2 2.97 1.48 13.3530 3.863e-06 ***
## G:T 2 0.10 0.05 0.4283 0.6522488
## M:T 2 0.02 0.01 0.0783 0.9247215
## A:T 2 0.06 0.03 0.2584 0.7725663
## F:G:M 1 0.05 0.05 0.4493 0.5035135
## F:G:A 1 0.03 0.03 0.3072 0.5800555
## F:M:A 1 0.36 0.36 3.2776 0.0718725 .
## G:M:A 1 0.03 0.03 0.2406 0.6243562
## F:G:T 2 0.82 0.41 3.6997 0.0266006 *
## F:M:T 2 0.05 0.02 0.2117 0.8093793
## G:M:T 2 0.05 0.03 0.2255 0.7983206
## F:A:T 2 0.09 0.04 0.3823 0.6828679
## G:A:T 2 0.54 0.27 2.4332 0.0905890 .
## M:A:T 2 0.63 0.32 2.8499 0.0604239 .
## F:G:M:A 1 0.01 0.01 0.0937 0.7598506
## F:G:M:T 2 0.09 0.04 0.3829 0.6824486
## F:G:A:T 2 0.05 0.02 0.2229 0.8004075
## F:M:A:T 2 0.57 0.29 2.5786 0.0786348 .
## G:M:A:T 2 0.12 0.06 0.5275 0.5909969
## F:G:M:A:T 2 0.10 0.05 0.4664 0.6280159
## Residuals 183 20.35 0.11
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(pH_anova)##
## Call:
## lm(formula = pH ~ F * G * M * A * T, data = abiotic)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.434 -0.102 -0.002 0.125 2.126
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.5600 0.1491 57.398 <2e-16 ***
## F -2.6420 0.2109 -12.527 <2e-16 ***
## G -0.0375 0.2237 -0.168 0.8671
## M 0.2660 0.2109 1.261 0.2088
## A 0.1860 0.2109 0.882 0.3790
## T24 0.0440 0.2109 0.209 0.8350
## T28 0.2480 0.2109 1.176 0.2412
## F:G 0.2315 0.3075 0.753 0.4524
## F:M -0.2920 0.2983 -0.979 0.3289
## G:M 0.0555 0.3075 0.181 0.8569
## F:A -0.1465 0.3075 -0.477 0.6343
## G:A -0.2065 0.3075 -0.672 0.5026
## M:A -0.7380 0.2983 -2.474 0.0143 *
## F:T24 0.1405 0.3075 0.457 0.6482
## F:T28 -0.0900 0.2983 -0.302 0.7632
## G:T24 0.2075 0.3075 0.675 0.5006
## G:T28 0.2375 0.3075 0.772 0.4408
## M:T24 -0.3900 0.2983 -1.308 0.1927
## M:T28 -0.1560 0.2983 -0.523 0.6016
## A:T24 -0.7500 0.2983 -2.515 0.0128 *
## A:T28 -0.3512 0.2983 -1.177 0.2405
## F:G:M 0.2425 0.4284 0.566 0.5720
## F:G:A 0.0330 0.4348 0.076 0.9396
## F:M:A 0.8865 0.4284 2.070 0.0399 *
## G:M:A -0.1315 0.4284 -0.307 0.7592
## F:G:T24 -0.3040 0.4348 -0.699 0.4853
## F:G:T28 -0.4455 0.4284 -1.040 0.2997
## F:M:T24 0.6215 0.4284 1.451 0.1485
## F:M:T28 0.0760 0.4218 0.180 0.8572
## G:M:T24 -0.2255 0.4453 -0.506 0.6132
## G:M:T28 -0.1135 0.4284 -0.265 0.7913
## F:A:T24 0.8305 0.4411 1.883 0.0613 .
## F:A:T28 0.3357 0.4284 0.784 0.4342
## G:A:T24 0.4165 0.4348 0.958 0.3394
## G:A:T28 0.2657 0.4284 0.620 0.5358
## M:A:T24 0.8760 0.4218 2.077 0.0392 *
## M:A:T28 0.5512 0.4218 1.307 0.1929
## F:G:M:A -0.3350 0.6058 -0.553 0.5809
## F:G:M:T24 -0.3505 0.6224 -0.563 0.5740
## F:G:M:T28 -0.1805 0.6012 -0.300 0.7643
## F:G:A:T24 -0.3390 0.6194 -0.547 0.5848
## F:G:A:T28 -0.2022 0.6058 -0.334 0.7389
## F:M:A:T24 -1.1965 0.6104 -1.960 0.0515 .
## F:M:A:T28 -0.6417 0.6012 -1.067 0.2872
## G:M:A:T24 0.2235 0.6179 0.362 0.7180
## G:M:A:T28 -0.1797 0.6012 -0.299 0.7653
## F:G:M:A:T24 0.4510 0.8749 0.515 0.6068
## F:G:M:A:T28 0.8202 0.8502 0.965 0.3360
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3335 on 183 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.9478, Adjusted R-squared: 0.9344
## F-statistic: 70.75 on 47 and 183 DF, p-value: < 2.2e-16oxygen_anova <- lm(Oxygen ~F*G*M*A*T, data = abiotic)
autoplot(oxygen_anova)
hist(residuals(oxygen_anova))
anova(oxygen_anova)## Analysis of Variance Table
##
## Response: Oxygen
## Df Sum Sq Mean Sq F value Pr(>F)
## F 1 24.59 24.594 6.9403 0.0091497 **
## G 1 21.31 21.315 6.0148 0.0151244 *
## M 1 47.88 47.878 13.5108 0.0003115 ***
## A 1 161.67 161.674 45.6227 1.845e-10 ***
## T 2 183.84 91.919 25.9385 1.209e-10 ***
## F:G 1 12.87 12.866 3.6306 0.0582953 .
## F:M 1 6.04 6.037 1.7036 0.1934602
## G:M 1 1.37 1.366 0.3856 0.5353968
## F:A 1 52.70 52.697 14.8705 0.0001593 ***
## G:A 1 4.33 4.334 1.2231 0.2702083
## M:A 1 0.23 0.231 0.0651 0.7989686
## F:T 2 34.02 17.012 4.8006 0.0092891 **
## G:T 2 9.56 4.779 1.3485 0.2622094
## M:T 2 3.51 1.756 0.4955 0.6100538
## A:T 2 50.91 25.454 7.1828 0.0009930 ***
## F:G:M 1 20.32 20.323 5.7349 0.0176401 *
## F:G:A 1 8.83 8.828 2.4913 0.1162046
## F:M:A 1 0.10 0.099 0.0278 0.8676644
## G:M:A 1 0.00 0.001 0.0003 0.9870868
## F:G:T 2 20.75 10.377 2.9283 0.0559984 .
## F:M:T 2 17.88 8.939 2.5225 0.0830454 .
## G:M:T 2 2.61 1.304 0.3681 0.6925778
## F:A:T 2 0.31 0.155 0.0437 0.9572609
## G:A:T 2 1.60 0.802 0.2264 0.7976449
## M:A:T 2 18.25 9.123 2.5745 0.0789470 .
## F:G:M:A 1 15.10 15.105 4.2624 0.0403759 *
## F:G:M:T 2 5.41 2.705 0.7634 0.4675410
## F:G:A:T 2 4.27 2.134 0.6021 0.5487344
## F:M:A:T 2 56.53 28.264 7.9758 0.0004774 ***
## G:M:A:T 2 2.23 1.114 0.3144 0.7306023
## F:G:M:A:T 2 14.74 7.371 2.0800 0.1278799
## Residuals 183 648.50 3.544
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(oxygen_anova)##
## Call:
## lm(formula = Oxygen ~ F * G * M * A * T, data = abiotic)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.2844 -0.7169 0.0606 0.7525 5.7695
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.3882 0.8419 23.030 < 2e-16 ***
## F 1.0142 1.1906 0.852 0.39541
## G 0.2158 1.2628 0.171 0.86450
## M -0.1926 1.1906 -0.162 0.87167
## A -2.6474 1.1906 -2.224 0.02740 *
## T24 -2.1344 1.1906 -1.793 0.07467 .
## T28 -1.1902 1.1906 -1.000 0.31879
## F:G -2.5408 1.7356 -1.464 0.14492
## F:M -1.6036 1.6837 -0.952 0.34215
## G:M -1.9018 1.7356 -1.096 0.27461
## F:A 3.2560 1.7356 1.876 0.06224 .
## G:A -1.1100 1.7356 -0.640 0.52326
## M:A 1.0388 1.6837 0.617 0.53803
## F:T24 4.5275 1.7356 2.609 0.00984 **
## F:T28 -1.0260 1.6837 -0.609 0.54304
## G:T24 1.9302 1.7356 1.112 0.26753
## G:T28 0.0004 1.7356 0.000 0.99982
## M:T24 1.4588 1.6837 0.866 0.38740
## M:T28 -0.2414 1.6837 -0.143 0.88615
## A:T24 -0.2418 1.6837 -0.144 0.88597
## A:T28 1.8400 1.6837 1.093 0.27591
## F:G:M 5.1164 2.4181 2.116 0.03571 *
## F:G:A 3.1290 2.4544 1.275 0.20399
## F:M:A -2.6114 2.4181 -1.080 0.28159
## G:M:A 3.6806 2.4181 1.522 0.12971
## F:G:T24 -5.3471 2.4544 -2.179 0.03064 *
## F:G:T28 0.6448 2.4181 0.267 0.79003
## F:M:T24 -4.8777 2.4181 -2.017 0.04514 *
## F:M:T28 0.7728 2.3812 0.325 0.74589
## G:M:T24 -2.0615 2.5139 -0.820 0.41325
## G:M:T28 2.1986 2.4181 0.909 0.36442
## F:A:T24 -5.9901 2.4903 -2.405 0.01715 *
## F:A:T28 -3.6768 2.4181 -1.521 0.13010
## G:A:T24 -1.4251 2.4544 -0.581 0.56221
## G:A:T28 0.9438 2.4181 0.390 0.69676
## M:A:T24 -2.8442 2.3812 -1.194 0.23384
## M:A:T28 -2.6162 2.3812 -1.099 0.27334
## F:G:M:A -6.1012 3.4197 -1.784 0.07606 .
## F:G:M:T24 2.8517 3.5134 0.812 0.41804
## F:G:M:T28 -4.6620 3.3937 -1.374 0.17121
## F:G:A:T24 3.3411 3.4965 0.956 0.34056
## F:G:A:T28 -1.5430 3.4197 -0.451 0.65237
## F:M:A:T24 9.2637 3.4455 2.689 0.00784 **
## F:M:A:T28 4.5152 3.3937 1.330 0.18502
## G:M:A:T24 0.7278 3.4881 0.209 0.83494
## G:M:A:T28 -5.4942 3.3937 -1.619 0.10718
## F:G:M:A:T24 -1.7131 4.9389 -0.347 0.72909
## F:G:M:A:T28 7.5692 4.7994 1.577 0.11650
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.882 on 183 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.5535, Adjusted R-squared: 0.4388
## F-statistic: 4.826 on 47 and 183 DF, p-value: 6.908e-15tc_anova <- lm(TC ~F*G*M*A*T, data = abiotic)
autoplot(oxygen_anova)
hist(residuals(oxygen_anova))
anova(oxygen_anova)## Analysis of Variance Table
##
## Response: Oxygen
## Df Sum Sq Mean Sq F value Pr(>F)
## F 1 24.59 24.594 6.9403 0.0091497 **
## G 1 21.31 21.315 6.0148 0.0151244 *
## M 1 47.88 47.878 13.5108 0.0003115 ***
## A 1 161.67 161.674 45.6227 1.845e-10 ***
## T 2 183.84 91.919 25.9385 1.209e-10 ***
## F:G 1 12.87 12.866 3.6306 0.0582953 .
## F:M 1 6.04 6.037 1.7036 0.1934602
## G:M 1 1.37 1.366 0.3856 0.5353968
## F:A 1 52.70 52.697 14.8705 0.0001593 ***
## G:A 1 4.33 4.334 1.2231 0.2702083
## M:A 1 0.23 0.231 0.0651 0.7989686
## F:T 2 34.02 17.012 4.8006 0.0092891 **
## G:T 2 9.56 4.779 1.3485 0.2622094
## M:T 2 3.51 1.756 0.4955 0.6100538
## A:T 2 50.91 25.454 7.1828 0.0009930 ***
## F:G:M 1 20.32 20.323 5.7349 0.0176401 *
## F:G:A 1 8.83 8.828 2.4913 0.1162046
## F:M:A 1 0.10 0.099 0.0278 0.8676644
## G:M:A 1 0.00 0.001 0.0003 0.9870868
## F:G:T 2 20.75 10.377 2.9283 0.0559984 .
## F:M:T 2 17.88 8.939 2.5225 0.0830454 .
## G:M:T 2 2.61 1.304 0.3681 0.6925778
## F:A:T 2 0.31 0.155 0.0437 0.9572609
## G:A:T 2 1.60 0.802 0.2264 0.7976449
## M:A:T 2 18.25 9.123 2.5745 0.0789470 .
## F:G:M:A 1 15.10 15.105 4.2624 0.0403759 *
## F:G:M:T 2 5.41 2.705 0.7634 0.4675410
## F:G:A:T 2 4.27 2.134 0.6021 0.5487344
## F:M:A:T 2 56.53 28.264 7.9758 0.0004774 ***
## G:M:A:T 2 2.23 1.114 0.3144 0.7306023
## F:G:M:A:T 2 14.74 7.371 2.0800 0.1278799
## Residuals 183 648.50 3.544
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(oxygen_anova)##
## Call:
## lm(formula = Oxygen ~ F * G * M * A * T, data = abiotic)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.2844 -0.7169 0.0606 0.7525 5.7695
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.3882 0.8419 23.030 < 2e-16 ***
## F 1.0142 1.1906 0.852 0.39541
## G 0.2158 1.2628 0.171 0.86450
## M -0.1926 1.1906 -0.162 0.87167
## A -2.6474 1.1906 -2.224 0.02740 *
## T24 -2.1344 1.1906 -1.793 0.07467 .
## T28 -1.1902 1.1906 -1.000 0.31879
## F:G -2.5408 1.7356 -1.464 0.14492
## F:M -1.6036 1.6837 -0.952 0.34215
## G:M -1.9018 1.7356 -1.096 0.27461
## F:A 3.2560 1.7356 1.876 0.06224 .
## G:A -1.1100 1.7356 -0.640 0.52326
## M:A 1.0388 1.6837 0.617 0.53803
## F:T24 4.5275 1.7356 2.609 0.00984 **
## F:T28 -1.0260 1.6837 -0.609 0.54304
## G:T24 1.9302 1.7356 1.112 0.26753
## G:T28 0.0004 1.7356 0.000 0.99982
## M:T24 1.4588 1.6837 0.866 0.38740
## M:T28 -0.2414 1.6837 -0.143 0.88615
## A:T24 -0.2418 1.6837 -0.144 0.88597
## A:T28 1.8400 1.6837 1.093 0.27591
## F:G:M 5.1164 2.4181 2.116 0.03571 *
## F:G:A 3.1290 2.4544 1.275 0.20399
## F:M:A -2.6114 2.4181 -1.080 0.28159
## G:M:A 3.6806 2.4181 1.522 0.12971
## F:G:T24 -5.3471 2.4544 -2.179 0.03064 *
## F:G:T28 0.6448 2.4181 0.267 0.79003
## F:M:T24 -4.8777 2.4181 -2.017 0.04514 *
## F:M:T28 0.7728 2.3812 0.325 0.74589
## G:M:T24 -2.0615 2.5139 -0.820 0.41325
## G:M:T28 2.1986 2.4181 0.909 0.36442
## F:A:T24 -5.9901 2.4903 -2.405 0.01715 *
## F:A:T28 -3.6768 2.4181 -1.521 0.13010
## G:A:T24 -1.4251 2.4544 -0.581 0.56221
## G:A:T28 0.9438 2.4181 0.390 0.69676
## M:A:T24 -2.8442 2.3812 -1.194 0.23384
## M:A:T28 -2.6162 2.3812 -1.099 0.27334
## F:G:M:A -6.1012 3.4197 -1.784 0.07606 .
## F:G:M:T24 2.8517 3.5134 0.812 0.41804
## F:G:M:T28 -4.6620 3.3937 -1.374 0.17121
## F:G:A:T24 3.3411 3.4965 0.956 0.34056
## F:G:A:T28 -1.5430 3.4197 -0.451 0.65237
## F:M:A:T24 9.2637 3.4455 2.689 0.00784 **
## F:M:A:T28 4.5152 3.3937 1.330 0.18502
## G:M:A:T24 0.7278 3.4881 0.209 0.83494
## G:M:A:T28 -5.4942 3.3937 -1.619 0.10718
## F:G:M:A:T24 -1.7131 4.9389 -0.347 0.72909
## F:G:M:A:T28 7.5692 4.7994 1.577 0.11650
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.882 on 183 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.5535, Adjusted R-squared: 0.4388
## F-statistic: 4.826 on 47 and 183 DF, p-value: 6.908e-15tn_anova <- lm(TN ~F*G*M*A*T, data = abiotic)
autoplot(tn_anova)
hist(residuals(tn_anova))
anova(tn_anova)## Analysis of Variance Table
##
## Response: TN
## Df Sum Sq Mean Sq F value Pr(>F)
## F 1 459470 459470 724.3140 <2e-16 ***
## G 1 12 12 0.0195 0.8891
## M 1 77 77 0.1212 0.7282
## A 1 1 1 0.0021 0.9634
## T 2 1964 982 1.5478 0.2155
## F:G 1 12 12 0.0194 0.8893
## F:M 1 132 132 0.2080 0.6489
## G:M 1 312 312 0.4919 0.4840
## F:A 1 173 173 0.2722 0.6025
## G:A 1 285 285 0.4493 0.5035
## M:A 1 143 143 0.2248 0.6360
## F:T 2 992 496 0.7818 0.4591
## G:T 2 129 65 0.1019 0.9032
## M:T 2 250 125 0.1970 0.8214
## A:T 2 2701 1351 2.1291 0.1219
## F:G:M 1 59 59 0.0932 0.7605
## F:G:A 1 54 54 0.0851 0.7708
## F:M:A 1 2 2 0.0038 0.9509
## G:M:A 1 56 56 0.0881 0.7669
## F:G:T 2 573 286 0.4513 0.6375
## F:M:T 2 139 70 0.1097 0.8962
## G:M:T 2 17 8 0.0132 0.9869
## F:A:T 2 1570 785 1.2376 0.2925
## G:A:T 2 223 111 0.1756 0.8391
## M:A:T 2 79 39 0.0622 0.9397
## F:G:M:A 1 8 8 0.0121 0.9124
## F:G:M:T 2 121 61 0.0957 0.9088
## F:G:A:T 2 176 88 0.1390 0.8704
## F:M:A:T 2 523 262 0.4126 0.6625
## G:M:A:T 2 195 98 0.1537 0.8576
## F:G:M:A:T 2 439 219 0.3457 0.7082
## Residuals 183 116086 634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(tn_anova)##
## Call:
## lm(formula = TN ~ F * G * M * A * T, data = abiotic)
##
## Residuals:
## Min 1Q Median 3Q Max
## -93.982 -0.366 0.085 1.879 61.348
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.21800 11.26368 0.019 0.985
## F 93.77400 15.92925 5.887 1.84e-08 ***
## G 0.31700 16.89552 0.019 0.985
## M 1.79800 15.92925 0.113 0.910
## A 0.08000 15.92925 0.005 0.996
## T24 0.11200 15.92925 0.007 0.994
## T28 0.20400 15.92925 0.013 0.990
## F:G -6.49300 23.22067 -0.280 0.780
## F:M -3.72000 22.52736 -0.165 0.869
## G:M 0.25100 23.22067 0.011 0.991
## F:A 0.16050 23.22067 0.007 0.994
## G:A -0.12100 23.22067 -0.005 0.996
## M:A -0.11800 22.52736 -0.005 0.996
## F:T24 -1.83400 23.22067 -0.079 0.937
## F:T28 -2.59600 22.52736 -0.115 0.908
## G:T24 -0.17100 23.22067 -0.007 0.994
## G:T28 0.24500 23.22067 0.011 0.992
## M:T24 0.10200 22.52736 0.005 0.996
## M:T28 -0.68800 22.52736 -0.031 0.976
## A:T24 16.93200 22.52736 0.752 0.453
## A:T28 -0.34800 22.52736 -0.015 0.988
## F:G:M 9.42700 32.35245 0.291 0.771
## F:G:A 7.01250 32.83898 0.214 0.831
## F:M:A 1.90950 32.35245 0.059 0.953
## G:M:A 0.29300 32.35245 0.009 0.993
## F:G:T24 10.56300 32.83898 0.322 0.748
## F:G:T28 8.94700 32.35245 0.277 0.782
## F:M:T24 -9.63000 32.35245 -0.298 0.766
## F:M:T28 5.18200 31.85849 0.163 0.871
## G:M:T24 -0.06033 33.63423 -0.002 0.999
## G:M:T28 -0.62700 32.35245 -0.019 0.985
## F:A:T24 -9.48250 33.31841 -0.285 0.776
## F:A:T28 2.31150 32.35245 0.071 0.943
## G:A:T24 -16.23950 32.83898 -0.495 0.622
## G:A:T28 -0.23900 32.35245 -0.007 0.994
## M:A:T24 -16.33400 31.85849 -0.513 0.609
## M:A:T28 0.78000 31.85849 0.024 0.980
## F:G:M:A -9.15450 45.75328 -0.200 0.842
## F:G:M:T24 -0.20617 47.00702 -0.004 0.997
## F:G:M:T28 -12.27500 45.40534 -0.270 0.787
## F:G:A:T24 8.34800 46.78157 0.178 0.859
## F:G:A:T28 -33.16650 45.75328 -0.725 0.469
## F:M:A:T24 26.51850 46.09859 0.575 0.566
## F:M:A:T28 -29.10350 45.40534 -0.641 0.522
## G:M:A:T24 16.12083 46.66843 0.345 0.730
## G:M:A:T28 0.00900 45.40534 0.000 1.000
## F:G:M:A:T24 -17.83783 66.07918 -0.270 0.788
## F:G:M:A:T28 35.37850 64.21284 0.551 0.582
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 25.19 on 183 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.8022, Adjusted R-squared: 0.7514
## F-statistic: 15.79 on 47 and 183 DF, p-value: < 2.2e-16#abiotic plots
abiotic <- abiotic %>%
reorder_levels(Combinations, order = c("Control","F","G","M","A","FG","FM","FA","GM","GA","MA","FGM","FGA","GMA","FMA","FGMA"))
stderr <- function(x, na.rm=FALSE) {
if (na.rm) x <- na.omit(x)
sqrt(var(x)/length(x))
}
#calculating mean of oxygen and pH
abiotic_means <- abiotic %>%
filter(T %in% c("20", "24", "28")) %>%
dplyr::group_by(Treatment, Combinations, T, Stressor,Stressor_numeric) %>%
dplyr::summarize(Oxygen_mean = mean(Oxygen),
pH_mean = mean(pH),
Oxygen_SE = stderr(Oxygen),
pH_SE = stderr(pH),
TN_mean = mean(TN),
TC_mean = mean(TC),
TN_SE = stderr(TN),
TC_SE = stderr(TC))## `summarise()` has grouped output by 'Treatment', 'Combinations', 'T', 'Stressor'. You can override using the `.groups` argument.plot_oxygen_means <- abiotic_means %>%
ggplot() +
geom_point(aes(x=T, y=Oxygen_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=T,y=Oxygen_mean,col=Combinations,ymin=Oxygen_mean-Oxygen_SE, ymax=Oxygen_mean+Oxygen_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
geom_line(aes(group=Combinations, x=T, y=Oxygen_mean, col=Combinations),position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
scale_colour_viridis_d(direction = -1)+
xlab("Temperature") +
theme(legend.position = "none") +
ylab("Oxygen[%] (mean)")
plot_oxygen_numeric_means <- abiotic_means %>%
ggplot() +
geom_point(aes(x=Stressor_numeric, y=Oxygen_mean, col=T, shape= T), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=Oxygen_mean,col=T,ymin=Oxygen_mean-Oxygen_SE, ymax=Oxygen_mean+Oxygen_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
geom_line(aes(group=Combinations, x=Stressor_numeric, y=Oxygen_mean, col=T),position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
ylab("Oxygen[%] (mean)") +
xlab("Number of stressors")+
theme(legend.position = "none")
#plot oxygen
plot_oxygen_means + plot_oxygen_numeric_means
#pH
plot_pH_means <- abiotic_means %>%
ggplot() +
geom_point(aes(x=T, y=pH_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=T,y=pH_mean,col=Combinations,ymin=pH_mean-pH_SE, ymax=pH_mean+pH_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
geom_line(aes(group=Combinations, x=T, y=pH_mean, col=Combinations),position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
scale_colour_viridis_d(direction = -1)+
xlab("Temperature") +
ylab("pH (mean)")+
theme(legend.position = "none")
plot_pH_means_numeric <- abiotic_means %>%
ggplot() +
geom_point(aes(x=Stressor_numeric, y=pH_mean, col=T, shape= T), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=pH_mean,col=T,ymin=pH_mean-pH_SE, ymax=pH_mean+pH_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
geom_line(aes(group=Combinations, x=Stressor_numeric, y=pH_mean, col=T),position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
ylab("pH (mean)")+
xlab("Number of stressors")+
theme(legend.position = "none")
#plot pH
plot_pH_means + plot_pH_means_numeric 
#total nitrogen
plot_TN_means <- abiotic_means %>%
ggplot() +
geom_point(aes(x=T, y=TN_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=T,y=TN_mean,col=Combinations,ymin=TN_mean-TN_SE, ymax=TN_mean+TN_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
geom_line(aes(group=Combinations, x=T, y=TN_mean, col=Combinations),position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
scale_colour_viridis_d(direction = -1)+
xlab("Temperature") +
ylab("TN (mean)")+
theme(legend.position = "none")
plot_TN_means_numeric <- abiotic_means %>%
ggplot() +
geom_point(aes(x=Stressor_numeric, y=TN_mean, col=T, shape= T), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=TN_mean,col=T,ymin=TN_mean-TN_SE, ymax=TN_mean+TN_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
geom_line(aes(group=Combinations, x=Stressor_numeric, y=TN_mean, col=T),position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
ylab("TN (mean)")+
xlab("Number of stressors")+
theme(legend.position = "none")
#plot nitrogen
plot_TN_means +theme(legend.position = "bottom",
legend.box = "vertical") + plot_TN_means_numeric + theme(legend.position = "bottom",
legend.box = "vertical")
#total carbon
plot_TC_means <- abiotic_means %>%
ggplot() +
geom_point(aes(x=T, y=TC_mean, col=Combinations, shape= Stressor), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=T,y=TC_mean,col=Combinations,ymin=TC_mean-TC_SE, ymax=TC_mean+TC_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
geom_line(aes(group=Combinations, x=T, y=TC_mean, col=Combinations),position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
scale_colour_viridis_d(direction = -1)+
xlab("Temperature") +
ylab("TC (mean)")+
theme(legend.position = "none")
plot_TC_means_numeric <- abiotic_means %>%
ggplot() +
geom_point(aes(x=Stressor_numeric, y=TC_mean, col=T, shape= T), size = 3, position = position_dodge(width = 0.4)) +
geom_errorbar(aes(x=Stressor_numeric,y=TC_mean,col=T,ymin=TC_mean-TC_SE, ymax=TC_mean+TC_SE), width=0.1,alpha=0.4,position = position_dodge(width = 0.4))+
geom_line(aes(group=Combinations, x=Stressor_numeric, y=TC_mean, col=T),position = position_dodge(width = 0.4)) +
theme_bw() + theme(panel.grid = element_blank()) +
theme(axis.text=element_text(size=12),
axis.title=element_text(size=14,face="bold"))+
ylab("TC (mean)")+
xlab("Number of stressors")+
theme(legend.position = "none")
#plot total carbon
plot_TC_means + plot_TC_means_numeric 
final plot for figure
(plot_diversity_means + plot_diversity_means_numeric)/
(nmds1_plot_temperature + plot_NMDS1_means_numeric)/
(nmds1_plot_temperature_subset + plot_NMDS1_means_numeric_subset)/
(plot_oxygen_means + plot_oxygen_numeric_means)/
(plot_pH_means +theme(legend.position = "bottom",
legend.box = "vertical",legend.text = element_text(size=12),legend.title = element_text(size=15)) +labs(col = "Stressor combinations",shape= "") + plot_pH_means_numeric + theme(legend.position = "bottom",
legend.box = "vertical",legend.text = element_text(size=12),legend.title = element_text(size=15))+ labs(col = "Temperature")+ labs(shape= "Temperature")) +plot_annotation(tag_levels = "a") 
#supplementary plot of all variables
(NMDS2_plot_temperature + plot_NMDS2_means_numeric) /
(NMDS2_plot_temperature_subset + plot_NMDS2_means_numeric_subset)/
(plot_TN_means + plot_TN_means_numeric) /
(plot_TC_means +theme(legend.position = "bottom",
legend.box = "vertical",legend.text = element_text(size=12),legend.title = element_text(size=15)) +labs(col = "Stressor combinations",shape= "") + plot_TC_means_numeric+ theme(legend.position = "bottom",
legend.box = "vertical",legend.text = element_text(size=12),legend.title = element_text(size=15))+ labs(col = "Temperature") + labs(shape= "Temperature")) + plot_annotation(tag_levels = "a")
extract the f and t test values, make one data frame for all factors
#NMDS1
ano_NMDS1_t <- tidy(anova_NMDS1)
names(ano_NMDS1_t)[names(ano_NMDS1_t) == "statistic"] <- "t_statistic"
names(ano_NMDS1_t)[names(ano_NMDS1_t) == "p.value"] <- "t_p.value"
names(ano_NMDS1_t)[names(ano_NMDS1_t) == "term"] <- "treatment_combination"
ano_NMDS1_t <- ano_NMDS1_t %>%
mutate(term = removeNumbers(treatment_combination))
ano_NMDS1_f <- tidy(anova(anova_NMDS1))
names(ano_NMDS1_f)[names(ano_NMDS1_f) == "statistic"] <- "f_statistic"
names(ano_NMDS1_f)[names(ano_NMDS1_f) == "p.value"] <- "f_p.value"
nmds1 <- full_join(ano_NMDS1_t, ano_NMDS1_f)## Joining, by = "term"nmds1 <- nmds1 %>%
mutate(color = ifelse(nmds1$f_p.value < 0.05, estimate, NA))
nmds1$variable="NMDS1"
nmds1 <- nmds1 %>%
filter(term != "(Intercept)", term!="Residuals")
nmds1 <- nmds1 %>%
reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))
#NMDS2
ano_NMDS2_t <- tidy(anova_NMDS2)
names(ano_NMDS2_t)[names(ano_NMDS2_t) == "statistic"] <- "t_statistic"
names(ano_NMDS2_t)[names(ano_NMDS2_t) == "p.value"] <- "t_p.value"
names(ano_NMDS2_t)[names(ano_NMDS2_t) == "term"] <- "treatment_combination"
ano_NMDS2_t <- ano_NMDS2_t %>%
mutate(term = removeNumbers(treatment_combination))
ano_NMDS2_f <- tidy(anova(anova_NMDS2))
names(ano_NMDS2_f)[names(ano_NMDS2_f) == "statistic"] <- "f_statistic"
names(ano_NMDS2_f)[names(ano_NMDS2_f) == "p.value"] <- "f_p.value"
nmds2 <- full_join(ano_NMDS2_t, ano_NMDS2_f)## Joining, by = "term"nmds2 <- nmds2 %>%
mutate(color = ifelse(nmds2$f_p.value < 0.05, estimate/1.9, NA))
nmds2$variable="NMDS2"
nmds2 <- nmds2 %>%
filter(term != "(Intercept)", term!="Residuals")
nmds2 <- nmds2 %>%
reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))
# NMDS1_CB_SB /1.8
ano_nmds1_mixed_t <- tidy(anova_NMDS1_mixed)
names(ano_nmds1_mixed_t )[names(ano_nmds1_mixed_t ) == "statistic"] <- "t_statistic"
names(ano_nmds1_mixed_t )[names(ano_nmds1_mixed_t ) == "p.value"] <- "t_p.value"
names(ano_nmds1_mixed_t )[names(ano_nmds1_mixed_t ) == "term"] <- "treatment_combination"
ano_nmds1_mixed_t <- ano_nmds1_mixed_t %>%
mutate(term = removeNumbers(treatment_combination))
ano_nmds1_mixed_f <- tidy(anova(anova_NMDS1_mixed))
names(ano_nmds1_mixed_f )[names(ano_nmds1_mixed_f ) == "statistic"] <- "f_statistic"
names(ano_nmds1_mixed_f )[names(ano_nmds1_mixed_f ) == "p.value"] <- "f_p.value"
nmds1_mixed <- full_join(ano_nmds1_mixed_t , ano_nmds1_mixed_f)## Joining, by = "term"nmds1_mixed <- nmds1_mixed %>%
mutate(color = ifelse(nmds1_mixed$f_p.value < 0.05, estimate/1.8, NA))
nmds1_mixed$variable="NMDS1_CB_SB"
nmds1_mixed <- nmds1_mixed %>%
filter(term != "(Intercept)", term!="Residuals")
nmds1_mixed <- nmds1_mixed %>%
reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))
#NMDS2_CB_SB
ano_nmds2_mixed_t <- tidy(anova_NMDS2_mixed)
names(ano_nmds2_mixed_t )[names(ano_nmds2_mixed_t ) == "statistic"] <- "t_statistic"
names(ano_nmds2_mixed_t )[names(ano_nmds2_mixed_t ) == "p.value"] <- "t_p.value"
names(ano_nmds2_mixed_t )[names(ano_nmds2_mixed_t ) == "term"] <- "treatment_combination"
ano_nmds2_mixed_t <- ano_nmds2_mixed_t %>%
mutate(term = removeNumbers(treatment_combination))
ano_nmds2_mixed_f <- tidy(anova(anova_NMDS2_mixed))
names(ano_nmds2_mixed_f )[names(ano_nmds2_mixed_f ) == "statistic"] <- "f_statistic"
names(ano_nmds2_mixed_f )[names(ano_nmds2_mixed_f ) == "p.value"] <- "f_p.value"
nmds2_mixed <- full_join(ano_nmds2_mixed_t , ano_nmds2_mixed_f)## Joining, by = "term"nmds2_mixed <- nmds2_mixed %>%
mutate(color = ifelse(nmds2_mixed$f_p.value < 0.05, estimate, NA))
nmds2_mixed$variable="NMDS2_CB_SB"
nmds2_mixed <- nmds2_mixed %>%
filter(term != "(Intercept)", term!="Residuals")
nmds2_mixed <- nmds2_mixed %>%
reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))
# Richness Shannon
ano_richness_t <- tidy(anova_diversity)
names(ano_richness_t)[names(ano_richness_t) == "statistic"] <- "t_statistic"
names(ano_richness_t)[names(ano_richness_t) == "p.value"] <- "t_p.value"
names(ano_richness_t)[names(ano_richness_t) == "term"] <- "treatment_combination"
ano_richness_t <-ano_richness_t %>%
mutate(term = removeNumbers(treatment_combination))
ano_richness_f <- tidy(anova(anova_diversity))
names(ano_richness_f )[names(ano_richness_f ) == "statistic"] <- "f_statistic"
names(ano_richness_f)[names(ano_richness_f ) == "p.value"] <- "f_p.value"
richness <- full_join(ano_richness_t, ano_richness_f)## Joining, by = "term"richness <- richness %>%
mutate(color = ifelse(richness$f_p.value < 0.05, estimate, NA))
richness$variable="Shannon_richness"
richness <- richness %>%
filter(term != "(Intercept)", term!="Residuals")
richness <- richness %>%
reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))
# Tychonema Cyanobacteria /2.2
ano_tychonea_t <- tidy(anova_tychonea)
names(ano_tychonea_t)[names(ano_tychonea_t) == "statistic"] <- "t_statistic"
names(ano_tychonea_t)[names(ano_tychonea_t) == "p.value"] <- "t_p.value"
names(ano_tychonea_t)[names(ano_tychonea_t) == "term"] <- "treatment_combination"
ano_tychonea_t <-ano_tychonea_t %>%
mutate(term = removeNumbers(treatment_combination))
ano_tychonea_f <- tidy(anova(anova_tychonea))
names(ano_tychonea_f )[names(ano_tychonea_f ) == "statistic"] <- "f_statistic"
names(ano_tychonea_f)[names(ano_tychonea_f ) == "p.value"] <- "f_p.value"
tyochonea <- full_join(ano_tychonea_t, ano_tychonea_f)## Joining, by = "term"tyochonea <- tyochonea %>%
mutate(color = ifelse(tyochonea$f_p.value < 0.05, estimate/2.2, NA))
tyochonea$variable="Tychonema"
tyochonea <- tyochonea %>%
filter(term != "(Intercept)", term!="Residuals")
tyochonea <- tyochonea %>%
reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))
# Sulfuricurvum /1.27
ano_sulfuri_t <- tidy(anova_sulfuri)
names(ano_sulfuri_t)[names(ano_sulfuri_t) == "statistic"] <- "t_statistic"
names(ano_sulfuri_t)[names(ano_sulfuri_t) == "p.value"] <- "t_p.value"
names(ano_sulfuri_t)[names(ano_sulfuri_t) == "term"] <- "treatment_combination"
ano_sulfuri_t <-ano_sulfuri_t %>%
mutate(term = removeNumbers(treatment_combination))
ano_sulfuri_f <- tidy(anova(anova_sulfuri))
names(ano_sulfuri_f )[names(ano_sulfuri_f ) == "statistic"] <- "f_statistic"
names(ano_sulfuri_f)[names(ano_sulfuri_f ) == "p.value"] <- "f_p.value"
sulfuri <- full_join(ano_sulfuri_t, ano_sulfuri_f)## Joining, by = "term"sulfuri <- sulfuri %>%
mutate(color = ifelse(sulfuri$f_p.value < 0.05, estimate/1.27, NA))
sulfuri$variable="Sulfuricurvum"
sulfuri <- sulfuri %>%
filter(term != "(Intercept)", term!="Residuals")
sulfuri <- sulfuri %>%
reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))
#oxygen /9.26
ano_oxygen_t <- tidy(oxygen_anova)
names(ano_oxygen_t )[names(ano_oxygen_t ) == "statistic"] <- "t_statistic"
names(ano_oxygen_t )[names(ano_oxygen_t ) == "p.value"] <- "t_p.value"
names(ano_oxygen_t )[names(ano_oxygen_t ) == "term"] <- "treatment_combination"
ano_oxygen_t <- ano_oxygen_t %>%
mutate(term = removeNumbers(treatment_combination))
ano_oxygen_f <- tidy(anova(oxygen_anova))
names(ano_oxygen_f)[names(ano_oxygen_f) == "statistic"] <- "f_statistic"
names(ano_oxygen_f)[names(ano_oxygen_f) == "p.value"] <- "f_p.value"
oxygen <- full_join(ano_oxygen_t , ano_oxygen_f)## Joining, by = "term"oxygen <- oxygen %>%
mutate(color = ifelse(oxygen$f_p.value < 0.05, estimate/9.26, NA))
oxygen$variable="oxygen"
oxygen <- oxygen %>%
filter(term != "(Intercept)", term!="Residuals")
oxygen <- oxygen %>%
reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))
#pH /2.62
ano_pH_t <- tidy(pH_anova)
names(ano_pH_t)[names(ano_pH_t) == "statistic"] <- "t_statistic"
names(ano_pH_t)[names(ano_pH_t) == "p.value"] <- "t_p.value"
names(ano_pH_t)[names(ano_pH_t) == "term"] <- "treatment_combination"
ano_pH_t <- ano_pH_t %>%
mutate(term = removeNumbers(treatment_combination))
ano_pH_f <- tidy(anova(pH_anova))
names(ano_pH_f )[names(ano_pH_f) == "statistic"] <- "f_statistic"
names(ano_pH_f )[names(ano_pH_f) == "p.value"] <- "f_p.value"
pH <- full_join(ano_pH_t , ano_pH_f)## Joining, by = "term"pH <- pH %>%
mutate(color = ifelse(pH$f_p.value < 0.05, estimate/2.62, NA))
pH$variable="pH"
pH <- pH %>%
filter(term != "(Intercept)", term!="Residuals")
pH <- pH %>%
reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))
#TN/93.77
ano_tn_t <- tidy(tn_anova)
names(ano_tn_t )[names(ano_tn_t ) == "statistic"] <- "t_statistic"
names(ano_tn_t )[names(ano_tn_t ) == "p.value"] <- "t_p.value"
names(ano_tn_t )[names(ano_tn_t ) == "term"] <- "treatment_combination"
ano_tn_t <- ano_tn_t %>%
mutate(term = removeNumbers(treatment_combination))
ano_tn_f <- tidy(anova(tn_anova))
names(ano_tn_f)[names(ano_tn_f) == "statistic"] <- "f_statistic"
names(ano_tn_f)[names(ano_tn_f) == "p.value"] <- "f_p.value"
tn <- full_join(ano_tn_t , ano_tn_f)## Joining, by = "term"tn <- tn %>%
mutate(color = ifelse(tn$f_p.value < 0.05, estimate/93.77, NA))
tn$variable="tn"
tn <- tn %>%
filter(term != "(Intercept)", term!="Residuals")
tn <- tn %>%
reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))
#TC /65
ano_tc_t <- tidy(tc_anova)
names(ano_tc_t )[names(ano_tc_t ) == "statistic"] <- "t_statistic"
names(ano_tc_t )[names(ano_tc_t ) == "p.value"] <- "t_p.value"
names(ano_tc_t )[names(ano_tc_t ) == "term"] <- "treatment_combination"
ano_tc_t <- ano_tc_t %>%
mutate(term = removeNumbers(treatment_combination))
ano_tc_f <- tidy(anova(tc_anova))
names(ano_tc_f)[names(ano_tc_f) == "statistic"] <- "f_statistic"
names(ano_tc_f)[names(ano_tc_f) == "p.value"] <- "f_p.value"
tc <- full_join(ano_tc_t , ano_tc_f)## Joining, by = "term"tc <- tc %>%
mutate(color = ifelse(tc$f_p.value < 0.05, estimate/65, NA))
tc$variable="tc"
tc <- tc %>%
filter(term != "(Intercept)", term!="Residuals")
tc <- tc %>%
reorder_levels(treatment_combination, order = c("T24","T28","F", "G","M","A","F:T24","G:T24", "M:T24","A:T24", "F:T28","G:T28","M:T28", "A:T28", "F:G", "F:M", "G:M", "F:A", "G:A", "M:A", "F:G:T24", "F:M:T24", "G:M:T24", "F:A:T24", "G:A:T24", "M:A:T24","F:G:T28", "F:M:T28", "G:M:T28", "F:A:T28", "G:A:T28", "M:A:T28", "F:G:M", "F:G:A","F:M:A","G:M:A","F:G:M:T24","F:G:A:T24","F:M:A:T24", "G:M:A:T24","F:G:M:T28", "F:G:A:T28","F:M:A:T28", "G:M:A:T28","F:G:M:A","F:G:M:A:T24", "F:G:M:A:T28"))
#all abiotic
all_abiotic <- rbind(oxygen, pH, tn, tc)
all_heat <- rbind(tyochonea,sulfuri,nmds1,nmds2,nmds1_mixed, nmds2_mixed,richness,oxygen, pH, tn, tc)
all_heat <- all_heat %>%
reorder_levels(variable,order=c("Tychonema","Sulfuricurvum","NMDS1","NMDS2","NMDS1_CB_SB","NMDS2_CB_SB","Shannon_richness","oxygen", "pH", "tn", "tc"))heat map of all factors
all_heat %>%
filter(term != "(Intercept)", term!="Residuals")%>%
ggplot(aes(x = variable, y = treatment_combination, fill=color)) +
geom_tile() +
scale_fill_gradient2(low="navy", mid="white", high="red",
midpoint=0,breaks=c(-1,0,1),
na.value="lightgrey",
limits=c(-1.2,1.2)
) +
theme(strip.placement = "inside") +
theme_bw()+
theme(axis.text.x = element_text(angle = 50, hjust = 1,size=10)) +
theme(axis.text.y = element_text(angle = 0, hjust = 1,size=10)) +
ylab("model term [main effect or interaction]")
just visualization of F:A with all factors
all_heat %>%
filter(treatment_combination %in% c("F","A","F:A","F:A:T24","F:A:T28"), variable != "tn",variable != "tc",color != "NA")%>%
ggplot()+
geom_point(aes(x=variable, y=color, shape=treatment_combination,col=treatment_combination),size=4)+
theme_bw() +
theme(axis.text=element_text(size=14),
axis.title=element_text(size=14)) +
ylab("standardised estimate")+
theme(axis.text.x = element_text(angle = 50, hjust = 1))+
geom_hline(yintercept=0)
numeric analysis
slope
#oxygen
abiotic<- ps_relative@sam_data %>%
as.tibble() %>%
select(Treatment,Temperature, Combinations, Temperature_factor, Lable, Name,F,G,M,A, Oxygen, TN,TC, pH,Stressor_numeric, Stressor)## Warning in class(x) <- c(setdiff(subclass, tibble_class), tibble_class): Setting
## class(x) to multiple strings ("tbl_df", "tbl", ...); result will no longer be an
## S4 objectabiotic <- abiotic %>%
filter(Treatment!="ori")
abiotic <- abiotic %>%
mutate(T2=Temperature^2,T=as.factor(Temperature)) %>%
select(-Temperature)
abiotic$Oxygen<- (abiotic$Oxygen-mean(abiotic$Oxygen))/sd(abiotic$Oxygen)
oxygen_anova_numeric <- lm(Oxygen ~T+T:Stressor_numeric, data=abiotic)
oxy_slope<- tidy(oxygen_anova_numeric ,conf.int = TRUE) %>%
filter(term!="(Intercept)", term!="T24", term !="T28")
oxy_slope$std.error <- NULL
oxy_slope$statistic <- NULL
oxy_slope$p.value <- NULL
oxy_slope$organisation="ecosystem"
oxy_slope$variable="oxygen"
#pH
abiotic$pH <- (abiotic$pH-mean(abiotic$pH))/sd(abiotic$pH)
pH_anova_numeric <- lm(pH ~T+T:Stressor_numeric, data=abiotic)
pH_slope <- tidy(pH_anova_numeric,conf.int = TRUE) %>%
filter(term!="(Intercept)", term!="T24", term !="T28")
pH_slope$std.error <- NULL
pH_slope$statistic <- NULL
pH_slope$p.value <- NULL
pH_slope$organisation="ecosystem"
pH_slope$variable="pH"
#TN
abiotic$TN <- (abiotic$TN-mean(abiotic$TN))/sd(abiotic$TN)
tn_anova_numeric <- lm(TN ~T+T:Stressor_numeric, data=abiotic)
tn_slope <- tidy(tn_anova_numeric,conf.int = TRUE) %>%
filter(term!="(Intercept)", term!="T24", term !="T28")
tn_slope$std.error <- NULL
tn_slope$statistic <- NULL
tn_slope$p.value <- NULL
tn_slope$organisation="ecosystem"
tn_slope$variable="total nitrogen"
#tc
abiotic$TC <- (abiotic$TC-mean(abiotic$TC))/sd(abiotic$TC)
tc_anova_numeric <- lm(TC ~T+T:Stressor_numeric, data=abiotic)
tc_slope <- tidy(tc_anova_numeric,conf.int = TRUE) %>%
filter(term!="(Intercept)", term!="T24", term !="T28")
tc_slope$std.error <- NULL
tc_slope$statistic <- NULL
tc_slope$p.value <- NULL
tc_slope$organisation="ecosystem"
tc_slope$variable="total carbon"
#NMDS1
mds_whole_res$NMDS1 <- (mds_whole_res$NMDS1-mean(mds_whole_res$NMDS1))/sd(mds_whole_res$NMDS1)
anova_NMDS1_numeric <- lm(NMDS1 ~T+T:Stressor_numeric, data=mds_whole_res)
nmds1_slope <- tidy(anova_NMDS1_numeric,conf.int = TRUE) %>%
filter(term!="(Intercept)", term!="T24", term !="T28")
nmds1_slope$std.error <- NULL
nmds1_slope$statistic <- NULL
nmds1_slope$p.value <- NULL
nmds1_slope$organisation="community"
nmds1_slope$variable="NMDS1"
#NMDS2
mds_whole_res$NMDS2 <- (mds_whole_res$NMDS2-mean(mds_whole_res$NMDS2))/sd(mds_whole_res$NMDS2)
anova_NMDS2_numeric <- lm(NMDS2 ~T+T:Stressor_numeric, data=mds_whole_res)
nmds2_slope <- tidy(anova_NMDS2_numeric,conf.int = TRUE) %>%
filter(term!="(Intercept)", term!="T24", term !="T28")
nmds2_slope$std.error <- NULL
nmds2_slope$statistic <- NULL
nmds2_slope$p.value <- NULL
nmds2_slope$organisation="community"
nmds2_slope$variable="NMDS2"
#NMDS1_subset
mds_whole_res_subset$NMDS1 <- (mds_whole_res_subset$NMDS1-mean(mds_whole_res_subset$NMDS1))/sd(mds_whole_res_subset$NMDS1)
anova_NMDS1_mixed_numeric <- lm(NMDS1 ~T+T:Stressor_numeric, data=mds_whole_res_subset)
nmds1_subset_slope <- tidy(anova_NMDS1_mixed_numeric,conf.int = TRUE) %>%
filter(term!="(Intercept)", term!="T24", term !="T28")
nmds1_subset_slope$std.error <- NULL
nmds1_subset_slope$statistic <- NULL
nmds1_subset_slope$p.value <- NULL
nmds1_subset_slope$organisation="community"
nmds1_subset_slope$variable="NMDS1_subpopulation"
#NMDS2_subset
mds_whole_res_subset$NMDS2 <- (mds_whole_res_subset$NMDS2-mean(mds_whole_res_subset$NMDS2))/sd(mds_whole_res_subset$NMDS2)
anova_NMDS2_mixed_numeric <- lm(NMDS2 ~T+T:Stressor_numeric, data=mds_whole_res_subset)
autoplot(anova_NMDS2_mixed_numeric)
nmds2_subset_slope <- tidy(anova_NMDS2_mixed_numeric,conf.int = TRUE) %>%
filter(term!="(Intercept)", term!="T24", term !="T28")
nmds2_subset_slope$std.error <- NULL
nmds2_subset_slope$statistic <- NULL
nmds2_subset_slope$p.value <- NULL
nmds2_subset_slope$organisation="community"
nmds2_subset_slope$variable="NMDS2_subpopulation"
#Shannon
div_raw$Shannon <- (div_raw$Shannon-mean(div_raw$Shannon))/sd(div_raw$Shannon)
anova_diversity_numeric <- lm(Shannon ~T+T:Stressor_numeric, data=div_raw)
autoplot(anova_diversity_numeric)
shannon_slope <- tidy(anova_diversity_numeric,conf.int = TRUE) %>%
filter(term!="(Intercept)", term!="T24", term !="T28")
shannon_slope$std.error <- NULL
shannon_slope$statistic <- NULL
shannon_slope$p.value <- NULL
shannon_slope$organisation="community"
shannon_slope$variable="species_richness"
#Tychonema
df_cyano$Abundance <- (df_cyano$Abundance-mean(df_cyano$Abundance))/sd(df_cyano$Abundance)
anova_tychonea_numeric <- lm(Abundance ~T+T:Stressor_numeric, data=df_cyano)
autoplot(anova_tychonea_numeric)
tychonema_slope <- tidy(anova_tychonea_numeric,conf.int = TRUE) %>%
filter(term!="(Intercept)", term!="T24", term !="T28")
tychonema_slope$std.error <- NULL
tychonema_slope$statistic <- NULL
tychonema_slope$p.value <- NULL
tychonema_slope$organisation="species"
tychonema_slope$variable="Tychonema"
#Sulfuri
df_sulfuri$Abundance <- (df_sulfuri$Abundance-mean(df_sulfuri$Abundance))/sd(df_sulfuri$Abundance)
anova_sulfuri_numeric <- lm(Abundance ~T+T:Stressor_numeric, data=df_sulfuri)
autoplot(anova_sulfuri_numeric)
sulfuri_slope <- tidy(anova_sulfuri_numeric,conf.int = TRUE) %>%
filter(term!="(Intercept)", term!="T24", term !="T28")
sulfuri_slope$std.error <- NULL
sulfuri_slope$statistic <- NULL
sulfuri_slope$p.value <- NULL
sulfuri_slope$organisation="species"
sulfuri_slope$variable="Sulfuricurvum"
all <- rbind(tychonema_slope,sulfuri_slope,shannon_slope,nmds1_slope,nmds2_slope,nmds1_subset_slope, nmds2_subset_slope,oxy_slope,pH_slope,tn_slope,tc_slope)
all <- all %>%
reorder_levels(variable,order=c("Tychonema","Sulfuricurvum","NMDS1","NMDS2","NMDS1_subpopulation","NMDS2_subpopulation","species_richness","oxygen","pH","total nitrogen","total carbon"))
plot_slope <-all %>%
ggplot(aes(x = organisation, y = estimate, colour = variable, shape=term)) +
geom_point(size = 4,position = position_dodge(width = 0.8)) +
geom_hline(aes(yintercept=0))+
geom_errorbar(aes(x=organisation,y=estimate,col=variable,ymin=conf.low, ymax=conf.high), width=0.1,alpha=0.4,position = position_dodge(width = 0.8))+
ylab("Slope [number of stressors]") +
theme_bw() +
scale_x_discrete(limits = c("species","community","ecosystem")) +
xlab("Level of organisation") +
theme(axis.title.x = element_text(size = 15), axis.text.x = element_text(size=12) ,axis.title.y = element_text(size = 15),axis.text.y = element_text(size=12))